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SYNOPSIS 


Automatic generation of the numerical grid is an important 
step in the computational simulation of a variety of applied 
problems. The requirements of automatic grid generation vary from 
problem to problem. In some cases, the grid is expected to 
smoothly fit a complicated geometry and yet have grid lines which 
intersect each other almost orthogonally (Thompson, 1984) . For 
situations which involve a moving front (say, a propagating flame 
or a shock wave) , capturing the location of the front and 
depicting the solution accurately near the front, become the 
primary objectives. In order to achieve these objectives, it is 
essential to generate the numerical grid in an adaptive fashion 
depending on the computed solution. 

In the present study, strategies for automatic grid 
generation have been developed based on the finite element 
approach, for both adaptive and non-adaptive situations. The 
applicability of these techniques has been illustrated on 
problems having complex geometries or moving fronts. 

The non-adaptive grid generation procedure developed in the 
present study is based on the solution of a Laplacian system of 
equations. It is meant for the generation of a structured, smooth 
grid inside an arbitrary shaped region. In this technique, the 
mapping of the complex physical domain into a simple shaped 
transformed domain is performed, with the help of body-fitting 
curvilinear coordinates (£,t?) (Fig. 1). 
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Fig 2 : Propagation of the flame front. 

A principle tool for such adaptive mesh generation is the 
development of an error estimator with a local a posteriori 
character, relating the accuracy of the computed solution to a 
perturbation of the original problem. An efficient a posteriori 
error estimate for finite element solutions has been derived in an 
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asymptotic form for step size h tending to zero, by Babuska and 
Rheinboldt (1978a) for a certain class of elliptic problems where 
the energy norm of the error is estimated as follows : 
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Here computable quantity tj is called the error indicator and is 
expressed by : 

,2 x 
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where n is the number of nodes in the domain, c is the constant 
and hj = x j+i - Xj. The a posteriori error in the equation (2) 
provides an easily computable bound for the total solution error 
under the energy norm. Also, the error indicator v provides a 
measure of error contribution of each element based on the local 
elemental residue r^ of the element j of the mesh. Using this 
concept of error analysis, Bieterman and Babuska (1982a) proposed 
an a posteriori error for the time dependent spatial meshes and 
presented the theory for the estimation of | |e(t,.) | | at each 
time. 


In the present study, the error estimates stated above are 
used in the h / hp-version adaption strategies. The a posteriori 
error estimate is defined by: 
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Here, Pj represents the local degree of finite element 

approximation considered in the j th element for the h- and 

hp-methods. The fundamental principle employed for grid adaption 

is the equi-distribution of the local error indicator i) within 

j 

the solution domain, while the a posteriori error £(t) is brought 
down below preset tolerance level 5 for each time. 

In the h-version method, the elements have been refined 
according to a statistical procedure based on the distribution of 
local error indicators (Oden, 1986b) . The level of refinement or 
de-refinement is selected on the basis of the deviation of the 
local error indicator from the mean value of the local error 
indicators. Such a statistical procedure involving simultaneous 
refinement and de-refinement at different parts of the solution 
domain is applied for each time step until the global a 
posteriori error meets the tolerance criterion. 

In the h-p method, the p-version and the h-version approaches 
are used alternately. For the p-strategy, the polynomial degree of 


interpolation in an 

element 

is 
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on 
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the 

deviation 
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mean of 

the 

error 
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the 
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is 


incremented (up to a maximum value of five ) and the global error 
criterion is checked for each adaption cycle. If the maximum 
allowable degree of interpolation is reached without satisfying 
the global error tolerance, then the concerned highest degree 
elements are divided according to the h-version strategy. In order 
to derive maximum benefit of the earlier p-version calculations, a 
gradient based non-uniform h-refinement has been introduced in 
the elements of highest degree approximation. 

Both the h- and hp-methods have been used to solve a test 
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problem with front propagation features and having a known 
analytical solution. Further, and the flame propagation problem 
given by equation (1) has also been solved. The exact and the 
computed solutions for the test problem have been compared 
extensively in respect of the solution profiles, the time-wise 
variations of quantities such as relative error, effectivity 
ratio, local and mean error indicators etc. , (Bieterman and 
Babuska, 1982a) . The h-version adaption method is seen to be more 
effective in bringing down the error than the hp-method. However, 
for the same global error tolerance, the number of elements 
required by the h-version is considerably higher. An elaborate 
comparison of the various error measures for the test problem 
indicates that the a posteriori error S(t) very closely represents 
the norm of the gradient error. 

Apart from the h- and hp-methods, a grid velocity scheme has 
been implemented in which the grid is moved with the speed of the 
front after the initial ignition phase. Such a grid speed strategy 
involves the least amount of computation for grid generation and 
it is also effective with respect to maintaining the solution 
error at a low value. But it is restricted to problems which have 
constant values of front width as well as front speed. 

A major contribution of the present work relates to the 
development of a generalized equi-distribution strategy 
(designated as the a-strategy) for the local error indicator and a 
posteriori error. Since the error estimators proposed by Bieterman 
and Babuska (1982a) are meant for self-adjoint spatial operators, 
their applicability to other types of problems, especially 
non-linear problems, is only tentative. The proposed a-strategy 
circumvents this difficulty. The only assumptions invoked for the 
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application of a- strategy are that the error estimator depends on 
the product of measure of the local residue in some properly 
chosen norm and positive power of the local mesh size. Based on 
these requirements the local error indicator used for the a 
posteriori error &(t) can be expressed as follows: 

2 a f Xj+1 2 

V ,(t) = h r. dx j = 0 , 1 , 2 , n-i (6) 

j J J x J 

j 

where the index a is itself an unknown and is to be determined 
after the calculation of the approximated solution. For any given 
mesh, the value of a is determined by equating the a posteriori 
error with a prescribed global tolerance limit 5. The degree of 
refinement for the j th element is decided according to the 
statistical procedure described earlier. Two limits of tolerance 
namely,: the upper and the lower tolerances and <5 2 ,are 

prescribed for the a posteriori error. The upper limit <5 1 is used 
to start the grid adaption, in the sense that whenever £(t) ^ 8^, 
grid is adapted. The lower limit is similarly used for stopping 
the adaption as well as for the determination of a. 

In order to illustrate the a-strategy, the problem of flame 
propagation in the presence of convective flow has been solved. 
This problem is described by the non-self -ad j oint operator of the 
form 


LI _ + y _i_T 

at ax 2 ax 
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= f(P.T) 
= - f(p.T) 


(7) 


where V is the fluid velocity. The boundary and initial conditions 
are prescribed as in equations (1-b) and (1-c) respectively. 

The Galerkin approach involving upwind weighting function has 
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been employed to obtain the finite element space discretization of 
the above equation, while the Crank-Nicolson scheme is employed 
for the purpose of time marching. The computed results have been 
compared with the analytical results of a test problem which 
displays the front movement feature. It is observed that the 
a-strategy performs exceedingly well in capturing the front 
movement and the spatial distribution of the local error indicator 
moves with front. Further, the computed error estimator approaches 
the norm of the actual gradient error, after a few initial cycles 
of grid adaption. For refined meshes, the value of a approaches 
two for a non-convective problem (V = 0) which agrees with the 
error estimate proposed by Bieterman and Babuska (1982a, b). The 
behavior of all the error measures indicates that the a-strategy 
is very effective in estimating and equi-distributing the solution 
error, in a self corrective fashion. 
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CHAPTER 1 

INTRODUCTION 


1.1. Motivation : 

Numerical solution of engineering problems has made great 
strides in recent years. Diverse fields such as aerodynamics, 
structural mechanics, process metallurgy etc. employ numerical 
simulation for analysing the optimality of a design or the 
efficiency of a process. A significant number of the application 
problems in these fields involve the solution of coupled partial 
differential equations in complex geometries. Over the past three 
decades, several numerical techniques have been developed for the 
solution of such partial differential equation (PDE) systems in 
regular domains. The focus of recent research is in the 
development of numerical techniques and software which facilitate 
the handling of complex geometries. A large volume of research 
literature has appeared on the generation of numerical grids in 
arbitrary- shaped geometries and various techniques have been 
proposed for generating meshes which are suitable for the 
application of finite difference (FD) , finite volume (FV) and 
finite element (FE) methods in solving the coupled PDE systems. A 
field of special interest is the dynamic generation of numerical 
grids for problems involving the propagation of sharp fronts. For 
instance, in premixed flame combustion, a flame front may 
propagate; in transient supersonic flows, moving shock fronts may 
occur. For the accurate capturing of such sharp fronts and for 
obtaining the overall solution to these problems, one needs to 
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take recourse to dynamic grid generation algorithms. 

For the development of grids in complex geometries, the use 
of body- fitting coordinate systems has found wide appea L in 
computational fluid dynamics and related fields. In order to 
resolve sharp fronts and the associated non-linear phenomena, 
solution- adaptive mesh generation techniques have been developed. 
During the generation of numerical grid by the adaptive or non- 
adaptive methods, several requirements are placed on the evolving 
grid structure. For instance, the grid lines should fit the 
geometry well at the boundary and vary smoothly within the 
interior of the domain. The grid should not become skewed and 
should not have too many or too few points at any particular 
region of the domain. Further, the grid spacing should be selected 
in such a way as to capture sharp fronts in the predicted 
solution. The grid spacing control is important not only for 
reducing the numerical error in the solution but also from the 
point of view of analysing important physical phenomena. Clearly, 
some of these requirements are mutually conflicting and one needs 
to find a compromise between them by judicial choice. Depending on 
the problem at hand and the numerical technique chosen to solve 
the PDE ' s , a particular grid characteristic is given more 
importance than the others. 

Although the developments of body- fitting grids for complex 
geometries and adaptive methods for problems with sharp fronts 
have advanced at a rapid pace in recent years, there are still 
several questions which are unanswered. Robust and quickly 
converging algorithms for the generation of smooth body- fitting 
grids with desired characteristics, are not yet available. 
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Further, the application of error- based adaptive strategies to 
non-linear problems such as flame propagation requires further 
study. In the present work, some finite element based non- 
adaptive and adaptive algorithms are developed and applied to 
arbitrary shaped geometries and flame propagation studies 
respectively. Detailed analysis of a posteriori error behaviour 
and the comparative features of alternative algorithms are 
evaluated. 

1.2. Survey of literature: 

Comprehensive surveys of literature related to the problems 
discussed in each of the subsequent chapters, have been presented 
in the respective chapters. Therefore, only an overview of the 
research carried out on adaptive and non- adaptive grid generation 
techniques is presented here. 

The grid generation techniques available at present can be 
classified into two major categories, namely: (i) Structured grid 
generation, ( i i ) Unstructured grid generation. Structured grid 
generation involves the transformation of a complex physical 
domain into a simple- shaped computational domain (usually 
rectangular) using curvilinear, body- fitting coordinates. The 
curvilinear coordinates are derived by one of the methods such as 
conformal mapping, algebraic interpolation or the solution of 
suitable elliptic partial differential equations. Structured grids 
are primarily useful for obtaining FD or FV based solutions. On 
the other hand, the unstructured grid generation techniques have 
been mainly used with finite element solution of the physical 
problem. These methods involve the algorithmic division of the 
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geometry into elements of chosen shape. Error analyses of the 
predicted FE solutions are often used during unstructured mesh 
generation for obtaining optimal mesh distributions in adaptive 
applications. Several excellent surveys, books and conference 
proceedings are available which summarize the research done on 
structured and unstructured grid generation ( Anderson, 1983b; 
Arcilla, 1991; Babuska and Rheinboldt, 1982; Babuska et al., 1983; 
Eiseman, 1986; Eiseman and Eriebacher, 1987; Ghia and Ghia, 1983; 
Hawken, 1987; Ho- Le, 1988 ; Sengupta et al . , 1988 ; Thompson et 
al., 1982, 1985a; Thompson, 1982a, 1984, 1985). A few important 

works on the different aspects of grid generation are described 
below. 

(a) Non- adaptive Methods for Structured Body- fitting Grids : 

In these methods, the complex geometry is transformed into a 
simple rectangular domain with the help of body- fitting 


coordinates 
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constant 
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movement of the physical boundaries. The transformation can be 
achieved analytically through conformal mapping ( O' Brien, 1981; 
Thompson et al., 1985a). However, conformal mapping is possible 
for regular geometries and two dimensional situations only. 
Numerical transformations can be obtained by algebraic 
(interpolation) techniques or through the solutions of elliptic 
partial differential equations. Algebraic methods are the fastest 
among all the grid generation techniques. They can be applied to 
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two- or three- dimensional problems ( Smith, 1983; Smith and 
Eriksson, 1980; Zhu et al . , 1988). However, boundary 
discontinuities are transmitted to the interior by these methods 
and therefore, obtaining smooth grids is difficult. The most 
general procedure for generating body- fitting grids involves the 
solution of elliptic PDE's for the variation of curvilinear 
coordinates ( Eiseman, 1985a; Ghia and Ghia, 1983; Mastin and 
Thompson, 1978; Thompson et al., 1977; Thompson, 1982b ). 
Smoothness and boundary orthogonality can be achieved in the grids 
generated by these methods. Indeed, with some compromise on 
smoothness and orthogonality, good control over grid spacing can 
also be obtained by selecting appropriate control functions in 
Poisson- type grid generation equations ( Steger and Sorenson, 
1979; Sorenson and Steger, 1980). The PDE grid solvers can be 
associated with variational functionals for grid smoothness, 
orthogonality and spacing control (Thompson et al., 1985a). The 
variational approach of obtaining structured grids has been 
discussed by Saltzmann and Brackbill (1982). 

(b) Non- adaptive Methods for Unstructured Grids : 

Unstructured meshes are generated by algorithmic division of 
the given geometry into elements of desired shape. These are 
associated with finite element solution methods. The available 
unstructured mesh generation techniques have been summarized by 
Thacker (1980) and Ho- Le (1988) . Generation of triangular 
elements have been discussed by Ho- Le (1988) and Joe & Simpson 
(1986) . Methods for obtaining quadrilateral elements in arbitrary 
2-D geometries have been described by Blacker and Stephenson 
(1991), Talbert and Parkinson (1990), Zhu et al . (1991) and 
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Zienkiewicz and Phillips (1971) . Quad-tree and Octree approaches 
have been developed which approximate the geometry in terms of 
step- like figures consisting of squares in 2-D and cubes in 3-D 
respectively (Kela et al., 1986; Yerry and Shepherd, 1984, 1988). 
Cavendish et al. (1984) have described another approach to three 
dimensional mesh generation . Geometrical methods of generating 
unstructured meshes have been developed by Bachmann et al. (1987), 
George et al. (1992) and Lo (1991). A versatile two dimensional 
mesh generator with automatic band width reduction has been 
proposed by Liu and Chen (1989). 

(c) Adaptive Techniques for Structured Grid Generation : 

It was described earlier that for elliptic grid solvers, 
spacing control is possible by proper choice of control functions 
(Sorenson and Steger, 1980; Steger and Sorenson, 1979; Brackbill, 
1982) . In adaptive methods, it is aimed to relate the spacing to 
some measure of local error and determine the spacing so that the 
spatial error is equi-distributed. The equi-distribution principle 
can be stated as 


p jCk, 

i+1 w(x) dx = constant 
J x 

i 

or, in discrete form, h t w = con stant, 

where h ][ (= x i + i ~ x t ) is the step size and w(x) is some positive 

weight function related to the magnitude of the local error in the 
predicted solution. Several adaptive grid methods have been 
devised which employ the equi- distribution concept and these are 
summarized in the review articles by Anderson (1983b), Berger 
(1987), Eiseman (1986), Ghia et al . (1983), Thompson (1985), 

Thompson et al. (1985b) . 
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Grid adaption based on the equi-distribution of solution 
gradients was studied by Dwyer et al. (1980, 1982), Dwyer and 
Sanders (1983) and Dwyer (1983). Nakahashi and Deiwert (1985) 
proposed a spring analogy approach to move the grid- points. The 
grid points were considered to be connected by hypothetical 
tension and torsion springs. By relating the spring constants to 
the solution gradients (or any other suitable measure of error) , 
the authors showed that the grid spacing and skewness could be 
effectively controlled. Saltzman and Brackbill (1982) and Kreis et 
al. (1986) have described the use of variational methods for 
adaptive grids. 

Movement of grid points with variable speeds until the 
equi-distribution principle is satisfied, has been suggested by 
many authors (Rai and Anderson, 1981; Anderson and Rai, 1982? 
Anderson, 1983a; Nakumara, 1982; Mack et al . , 1992). Rai and 
Anderson (1981) developed the attraction- repulsion strategy for 
evaluating the grid speed. The nodal speeds were determined on the 
basis of local truncation error. Nakumara (1982) proposed the use 
of parabolic partial differential equations for moving the grid 
points . 

(d) Adaptive Methods for Unstructured Grids : 

The research literature in this area has been summarized by 
Babuska and Rheinboldt (1982), Babuska et al. (1983), Carey and 
Humphrey (1981) , Carey and Oden (1984a) , Thompson (1985) and Zhu 
and Zienkiewicz (1991) . 

Adaptive triangulation procedures have been described by 
Eiseman (1985b) , Eriebacher and Eiseman (1987) and Rivara (1984) . 
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Shepherd et al. (1986) and Yerry and Shepherd (1984) have 
developed octree based adaptive mesh refinement. Hierarchical FE 
refinement has been discussed by Adjedj and Aubey (1989) and Yuen 
et al. (1991). Moving finite elements have been studied by Adjerid 
and Flaherty (1986a, b), Arney et al . (1986), Gelinas et al . (1981) 
and Miller and Miller (1981a, b). Adaptive time discretization has 
been considered by Babuska and Luskin (1989) and the space- time 
approach has been implemented by Bjer (1989) . 

The most significant development in adaptive FE procedures is 
the construction of a posteriori estimates for the spatial error 
in the computed solution and grid adaption based on the equi- 
distribution of local errors. The research in this area has been 
pioneered by Babuska and Rheinboldt (1978a-b, 1979a-c, 1980, 1981, 
1982), Bank & Weiser (1985), Bank (1986), Bieterman and Babuska 
(1982a-b, 1986), Carey (1988), Demkowicz et al. (1989), Flaherty 
et al. (1983) ,Gago et al. (1983) and Kelly et al (1983). Based on 
these estimates, h- version methods which refine the grid spacing 
and p- version techniques which increment the order of FE 
interpolation have been developed by Babuska and Dorr (1981), 
Babuska et al. (1981), Babuska et al. (1984), Demkowicz et al . 
(1989), Gui and Babuska (1986a-c) and Rachowicz et al. (1989). The 
concept of feedback for grid adaption has been introduced by 
Rheinboldt (1983), Babuska (1986) and Babuska and Gui (1986). 

( e ) Applications of Adaptive Grids : 

The application which is of primary interest to the present 
study is the numerical modelling of combustion and the associated 
flame propagation phenomena. Several studies have been carried out 
in recent years on this subject. These have been pioneered by 
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Benkhaldoun and Larrouturou (1987, 1989), Benkhaldoun et al. 

(1988a-b) , Dervieux et al. (1989), Dervieux and Larrouturou 
(1990), Dwyer and Sanders (1983), Glowinski (1985), Larrouturou 
(1985, 1986) and Ramos (1983, 1985, 1987, 1990a-b, 1991). Various 

aspects of combustion and flame propagation have been analysed by 
the above authors. A majority of these studies employ finite 
difference based grid adaption or the moving mesh methods. 
Incorporation of h-version FE adaption has also been carried out 
by some authors. 

1.3 Objectives of the Present Work : 

One of the objectives of the present study is to simplify the 
solution procedure for the elliptic equations used in the 
generation of structured grids in complex geometries. Although the 
elliptic PDE based grid generation systems are versatile, their 
applicability at present is being limited by the strong non- 
linearity of the transformed equations in the computational 
domain. Here, it is proposed to solve the elliptic PDE's in the 
physical domain itself using the finite element method, without 
transforming them into the rectangular computational domain. By 
such simplification, a robust and versatile grid generation 
algorithm can be constructed for obtaining structured grids with 
the desired level of smoothness, orthogonality and spacing 
control . 

Another major objective is to model one- dimensional flame 
propagation using the error estimates derived by Babuska and 
Rheinboldt (1982) in the h- and h-p version refinement algorithms. 
It is intended to compare the relative performances of the two 
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methods among themselves and also against the moving mesh methods, 
for the flame propagation problem . A detailed study of the error 
behaviour and the optimality of grid distribution, will also be 
carried out. 

The final objective is to develop a generalized self- 
correcting error estimator strategy which eliminates the need for 
a detailed mathematical analysis for constructing appropriate a 
posteriori error estimates. It is intended to test this algorithm 
on the flame propagation problem which includes fluid flow. 

1. 4 Outline of the Dissertation : 

In the Chapter 2, structured grid generation using the finite 
element method is described. The method has been applied for the 
generation of body- fitting grids in simply and multiply connected 
domains. The effects of factors such as the nature of boundary 
data specified and the initial guess grid, on the characteristics 
of the converged grid are examined. 

In Chapter 3, the one- dimensional flame propagation during 
premixed combustion is described. The mathematical framework for 
the grid adaption study and the definition of a posteriori error 
estimates and associated quantities, are provided. The finite 
element solution procedure for the flame propagation in a 
stationary medium is discussed. 

In Chapter 4, the h- , h-p, and grid speed strategies are 
described. All the steps involved in the grid refinement or the 
derefinement procedures are elaborately discussed. The predicted 
results for the three adaption strategies are compared for both 
the test and the flame problem. The computed results are also 
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validated with the help of a test problem which has similar front 
movement features. 

In Chapter 5, the generalized self- correcting error 
estimator strategy is described. Its application is illustrated 
with the help of the test problem and the flame problem, for a 
flowing mixture. 

In Chapter 6, the main conclusions of the work are summarized 
and some suggestions for future work are provided. 



CHAPTER 2 


GENERATION OF BODY-FITTING GRIDS USING FINITE ELEMENT METHOD 

2.1. Introduction : 

In this chapter, a novel procedure for the automatic 
numerical generation of a structured, smooth grid with 
body-fitting coordinates inside an arbitrarily shaped region is 
discussed. The method is based on Laplace grid generation 
equations. The finite element solution procedure and the 
construction of body conforming coordinates (£, tj) , have been 
discussed. The method has been applied for generating grids in 
simply and multiply connected two dimensional regions. Detailed 
numerical experimentation has been performed on the construction 
of O-type grids, and the effects of various features such as the 
density of boundary grid point data, the shape of the branch cut 
and the starting guess grid, have been highlighted. 

2.2. Literature survey: 

Numerical solutions of most engineering problems involve the 
use of Finite Difference (FD) or Finite Element (FE) methods. For 
both the methods, the first step concerns with the discretization 
of the spatial domain into a suitable numerical grid. There are, 
however, differences in the required grid structure for the two 
methods . 

For FD applications, it is most convenient to consider a 
regularly placed, rectangular array of grid points. Therefore, a 
problem with complex geometry is tackled by mapping the 
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physicaldomain into a rectangular computational domain, by 
suitable transformation. Body fitting co-ordinates are employed 
for transforming the geometry and the implementation of boundary 
conditions becomes simple when expressed in terms of the 
body-fitting co-ordinates. Additional requirements on the grid 
structi^ for these applications are that the spacing between the 
grid lines must be controllable (based on the gradients of the 
solution variables) and the grid lines should not become too 
skewed. 

For FE methods, the role of automatic grid generation is to 
simplify the specification of input data to the numerical 
simulation software (pre-processing) . It is necessary to number 
the nodes and elements of the mesh and also generate the data of 
nodal co-ordinates and the nodal connectivity of elements. Here, 
the complex geometry of a given solution domain is divided 
algorithmically into elements of required shape and the other grid 
data are automatically computed by some algorithms. The grid 
generation techniques available at present fall into two 
categories, namely: a) structured grid generation b) unstructured 
grid generation. Structured grid generation produces a regular 
array of points and it is suitable for the transformation of a 
complex geometry into a rectangular one by the use of body-fitting 
co-ordinates. Therefore, structured grid generation is widely used 
along with FD simulation software. Unstructured grid generation on 
the other hand, has no constraints such as the production of a 
regular array of points. The basic aim here is to approximate a 
complex geometry in terms of elements of chosen shape, in the best 
manner possible. The unstructured grid generation, is suited for 
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FE applications. The earlier work on these generation schemes have 
been briefly discussed below. 

2 . 2 . 1 . Unstructured grid generation : 

Most of the work in unstructured mesh generation relates to 
techniques such as 2-D quadree and 3-D Octree approaches, and 
triangulation or sub structuring of the geometry ( Batchman et 
al., 1987 ; Kela et al., 1986 ; Shepherd and Yerry, 1983; Joe and 
Simpson, 1986 ; Cook, 1981 ). One of the widely adopted techniques 
for unstructured grid generation is the automatic triangulation, 
which involves algorithmic division of the physical region into 
nearly equilateral triangular elements. Some of these techniques 
have been extended to produce all-quadrilateral meshes ( Shepherd, 
1984 ; Thacker, 1980 ). Several authors employ parameter space 
mapping to generate quadrilateral meshes interactively (Cook, 
1981 ; Gordan; 1973 ; Ho-Le, 1988 ). Techniques for the local 

adaptive refinement of the element sizes in required sub regions 
have also been developed ( Yerry and Shepherd, 1988; Zhu and 
Zienkiewicz, 1991 ) . 

2.2.2. Structured grid generation : 

The structured grid generation techniques are based on the 
transformation of the complex physical domain into a simple 
computational domain, which is often chosen to be rectangular in 
shape. Body-fitted curvilinear co-ordinates ( say £ and 7 }, where 
both £ and t; are functions of the physical co-ordinates x and y) 
are used for the transformation. The numerical grid is developed 
by plotting the curvilinear co-ordinate lines i.e., £ = constant 
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and 7? = constant curves. The governing equations of the problem to 
be solved are converted in terms of the curvilinear co-ordinates 
and the resulting system of equations are solved in the 
transformed (rectangular) domain by finite difference methods 
(FDM) . The transformation into the domain results in an 
orthogonal and uniform computational grid which is required by 
FDM. Its advantage lies in the relative easiness of programming 
since all the calculations are on the regular grid. Moreover, 
curved boundaries are accounted for and there is no loss of 
accuracy due to inaccurate modeling of boundary. Further, the grid 
point positions can be easily changed in the physical domain 
according to the needs of the computed solution. 

Two important questions arise during the generation of body 
conforming grids. First, how to automatically generate grids which 
are suited for the problem at hand, and, second, how to control 
the resulting grid point distribution. In order to address these 
questions a variety of techniques have been used from simple 
co-ordinate stretching to the solution of elliptic equations for 
the curvilinear co-ordinates ( Eiseman and Erlebacher; 1987 
Thompson, 1982b, 1984, 1985a-d ; Steger and Sorenson, 1979; Smith 
and Eriksson, 1986 ) . Among these, three important techniques 
which deserve attention are the algebraic grid generation, 
conformal mapping and the elliptic PDE solvers. For open 
(external) domains, parabolic and hyperbolic grid generation 
methods have been developed using marching techniques ( Nakumara, 
1982; Dwyer et al., 1982 ). 

A majority of algebraic grid generation techniques employ 
fast interpolation methods. Among the various interpolation 
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procedures, transfinite interpolation is a widely used method. In 
transfinite interpolation, the transformation procedure involves 
the use of a curvilinear coordinate system which conforms to the 
shape of the physical region at the boundaries. The interior 
coordinate lines are obtained using interpolation between the 
boundary values of the co-ordinates. Grid spacing can also be 
controlled by prescribing the step size variation on fixed 
boundaries according to some algebraic expressions (Eiseman, 1986; 
Eiseman and Erlebacher; 1987; Smith and Erikkson; 1986) . These 
methods are very simple to use. However, for complex geometries, 
it is difficult to represent the boundary surfaces accurately and 
the interpolation technigues fail if some parts of the boundary 
are described by concave surfaces. Handling of sharp corners on 
the boundary is also not satisfactory. 

Complex variable methods which involve conformal mapping have 
also been widely used (Eiseman and Erliebacher, 1987) . 

Unfortunately they are restricted to two dimensions only and also 
arbitrary shaped geometries cannot be handled. These techniques 
have been elaborately discussed in the survey papers, proceedings 
and the books published by various authors ( Eiseman, 1986; 
Thompson, 1982a, b, 1984; Thompson et al., 1982, 1985a-c; Turkel, 

1983; Hawken, 1987). 

Thompson et al., ( 1977, 1982, I985a-c ) have discussed the 
use of partial differential equations to generate grids in 
arbitrary geometries. The PDE grid solvers can meet most of the 
important requirements such as smoothness and orthogonality of 
grid lines, and control of grid-spacing. Laplace grid generation 
equations have been .considered for obtaining the body-fitting 
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coordinates during the geneation of smooth grids (Thompson et. 
al., 1985a). These equations have been solved by the finite 
difference methods over the transformed domain (rectangular) 
with physical variables x and y as the unknowns. The above 
authors have derived the orthogonal families of body-fitting 
curves by solving the Laplace equations. The orthogonal 
body-fitting grid gives a better representation of a general 
complex geometry because of its ability to fit a grid line to each 
boundary surface. However, the transformed equations become quite 
complex and there is only limited control over the refinement of 
the generated grid because of the need to retain orthogonality. 
If orthogonality is sacrificed, then body fitting grids may be 
derived from the solutions of the Poisson systems. The right hand 
side of the Poisson equations exercise control over the 
clustering of grid lines or grid points which is very important 
from the applications point of view. The use of such control 
functions on the right hand side allows a very fine resolution at 
certain specified locations and it covers adaptive grid generation 
as well. On the other hand, there is a chance that the extremum 
principles of the Laplace operator may be lost due to the choice 
of control functions. Many choices for the construction of 
control functions are discussed in the works of Thompson et al., 

( 1977, 1982, 1985a-c) and Eiseman (1986). 

During the transformation of the physical partial 
differential equations, the gradient, divergence, curl, and 
Laplace operators get replaced by either their conservative or 
non-conservative expressions in terms of general curvilinear 
co-ordinates. Thompson et al., (1985a) have presented the 
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transformed derivatives ( metric coefficients) using the notion of 
co-variant and contravariant directional derivatives and vector 
components. From the truncation error point of view, the metric 
coefficients have to be evaluated accurately. This is because they 
appear as coefficients for the transformed governing equations of 
the physical problem as well as the Poisson grid generation 
system. Thus, grid generation using elliptic PDE is a complex, 
but very general procedure for obtaining body-fitting grids, with 
appropriate control for achieving desired grid characteristics. 
Since the present study deals with the generation of body-fitting 
grids using an elliptic PDE system, we briefly describe this 
procedure in the following section. 

2.2.3. Grid generation equations based on elliptic equations : 

In this subsection, the general two dimensional coordinate 
transformation procedure and the transformed grid generation 
equations are presented (Thompson et al. 1985a). These have been 
widely used in the past for finite difference simulation of 
physical problems with complex geometries. 

Consider a simply or multiply connected region Q(x,y) £ (R 2 of 
arbitrary shape to be transformed into a rectangular region 
as shown in the Figs. 2.1a-b. The general transformation 

from the physical plane (x,y) to the transformed plane (£, 17 ) is 
given by 

? = £(x,y), 7 ) = 7 ] (x, y) . ( 2 . 1 a) 

Similarly, the inverse transformation is given by 

x = x(£,t 7) , y = y(£,v) . (2.1b) 

For a given function f = f (x,y) , the transformed derivatives can 


A 
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Fig. 2.1 (a) Transformation of an ellipse to a rectangle. 



Fig. 2.1(b) Transformation of an annular region to a rectangle. 
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be written as follows : 


d(f,y ) d(x,y) 

/ 


d(Z,V) 


d(Z'V) 


d(fi'Tl) 


ay af ay ('){' 
dr, 3£ 3£ dr) 


d(f,x) d(x,y) 

/ 


3(€,t}) 




/ J 


dx af ax af . . 
dr) 3£ 3£ ar) I 7 


( 2 - 2 ) 


j _ dx 3y 


dx ay 
dr) 377 


where J is the Jacobian of transformation given by , 

(2.3) 

The basic idea of the transformation is to generate transf ormation 
functions such that all boundaries are coincident with co-ordinate 
lines. The body-fitting co-ordinates (£ , 77 ) are taken as the 
solutions of a suitably chosen elliptic boundary value problem, 
with one of the co-ordinates constant on the boundaries. The 
Poisson equations associated with grid generation can be 
written as 


+ Ik 

a 2 x' a 2 y > 

8% + 3J1 

a 2 x x a 2 y 


P(C/T7) 

Q(ZrV) 


(2.4a) 

(2.4b) 


where x - x{£,r}) , y = y(£,77) and x,y e Q e CR 2 with suitable 
boundary conditions which specify either the grid spacing or the 
grid slope at the boundary. Interchanging dependent and 
independent variables for eqns. (2.4a-b) gives : 


a 

ae 2 

- 2/3 

a 2 x + 

d^dr, 

,2 

r 8 J 

377 

+ j 2 (p || 

+ Q ) - 0 

(2.5a) 

a 2 y 
a — 

a£ 2 

- 2/3 

d 2 y + 

df[dT) 

a 2 y 

* 2 

dr, 2 

+ j 2 ( p §| 

+ Q |f , - 0 

(2.5b) 
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where a = 






3x Sx ay ay 
a£ 5r) + a£ an nd 






The equations (2.5a-b) can now be solved on the transformed £- 7 7 
plane using finite difference methods, and the values of the 
physical co-ordinates (x, y) for each grid point (£, tj) in the 
transformed plane can be obtained. It is seen that the equations 
in the transformed plane are more complicated than in the physical 
plane, however with the advantage that the solutions are being 
found for a simple square grid. Use of the boundary-fitted 
co-ordinate system enables the boundary conditions to be applied 
easily. 


2.3. Scope of present work : 

In the previous section, the generalized approach of elliptic 
PDE based grid generation procedure was outlined. It can be shown 
(Thompson et. al. 1985; Saltzman and Brackbill, 1982) that the 
Poisson grid generation equations correspond to the minimization 
of a functional which is made up of terms related to grid 
smoothness, orthogonality and spacing control. In this sense, the 
elliptic PDE solvers can provide optimal grids, for a given 
physical problem. However, as discussed in the above section, when 
the grid generation equations are transformed in terms of the body 
fitting coordinates for facilitating FD solutions, the resulting 
equations are highly nonlinear and complex. These have to be 
solved iteratively and the iterative procedure often poses serious 
difficulties for convergence. 

In the present work, structured grid generation by the finite 



element solution of elliptic partial differential equations has 
been studied! The generation of structured grids in the physical 
domain by finite element method can offer a lot of advantages. 
This is due to fact that the original Poisson grid generation 
equations (2.4 a, b) are linear and these can be solved easily in 
the physical domain itself . Consequently important grid features 
such as orthogonality, smoothness etc. can be preserved, since no 
non-linear iterations are involved. Most numerical algorithms are 
very sensitive to grid skewness, clustering and smoothness; if 
grid generation procedure does not converge properly, the 
evaluation of metric terms may affect the numerical accuracy of 
the predicted solutions for such schemes. Thus, there is a strong 
need to simplify the solution procedure for the elliptic PDE 
solvers and develop rapidly convergent grid generation schemes. 

The generation of smooth body fitting grids in complex 

geometries using Laplace equations is considered in the present 
* 

work. Since all the computations are performed in the physical 
domain only, there is no need to transform the grid generation 
equations in terms of the curvilinear co-ordinates. Apart from 
linearity, the other inherent advantages of the Laplacian grid 
solver are also intactly preserved. The implementation of 
Dirichlet and Neumann boundary conditions which relate to the 
grid spacing and the slope of grid lines at the boundary 
respectively, is easily achieved in the proposed scheme. 

2.4. Geometry of the problems considered : 

Description of the Geometry : 

The present method has been applied to a simply connected 

This work is published in the conference (see, Rao et al, 1988). 
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domain (Fig. 2.1a) and a multiply connected region with an O-type 
grid (Fig. 2.1b). In the former case, the transformation is from 
an ellipse to a rectangular region and in the latter, it is from 
an annular region to a rectangle with the use of a branch cut. 
Different grid configurations of the body-fitting co-ordinate 
system can be created by identifying different boundary segments 
of the physical region with paritcular sides of the transformed 
rectangular region( Thompson et al., 1985a). 

Boundary conditions : 

The boundary conditions for the elliptic domain problem are 

of Dirichlet type (Fig. 2.1a). In this case the boundary of the 

ellipse is divided into two £- constant curves ( C 2 and C 4 ) on 

which T) varies monotonically and two ij- constant curves ( C and 

C ) on which £ varies monotonically. For the annular region (Fig. 

2.1b), the £- variation is specified on the i n constant curves C 

and C 3 . On the £- constant curves C 2 and C^, the 7)-variation is 

left unspecified and re-entrant boundary conditions for r) are 

enforced. In the finite element implementation, the re-entrant 

boundary conditions are applied very easily by declaring the 

respective elements and nodes across the branch cut to be 

neighbours (by specifying the nodal connectivity of elements 

appropriately) . Apart from Dirichlet boundary conditions for £ on 

C i and C , Neumann boundary conditions (the gradient of £ in the 

direction normal to the boundary, |^) can also be specified 

depending upon the need. A particular case of interest is the 

orthogonal grid at the boundary; for this, the value of the 

gradient, 4^ is taken as zero. 
on 
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2.5. Finite element formulation of grid generation equations : 

The grid generation equations which have been used in the 
present work are of the form. 


a 2 £ 

a x 2 


_a%_ 

a x 2 


+ 


a 2 g 

a y 2 

j9JL 
a y 2 


= 0 


(2.6a) 


(2.6b) 


where x = x(£,7i) , y = y(£,T7) and x,y e Q £ IR , with the boundary 
conditions specified as mentioned above. 


Finite element discretization : 

The given physical domain Q(x,y) is divided into several 
small quadrilateral elements (Figs 2.2 and 2.3). Eight-noded 
curvilinear elements have been employed in the present work, which 
can fit arbitrary geometries quite well. Isoparametric mapping of 
each curvilinear element has been performed with the help of local 
curvilinear coordinates x and y, as shown in Fig. 2.2. As per 
conventional FE procedure, such local mapping facilitates 
numerical quadrature during the calculation of elemental matrices, 
corresponding to the discretized forms of the grid generation 
equations (2.6a,b). 

The construction of the numerical grid using £- and 
7)-constant curves can be explained in the following manner. We 
assume that on any £ curve for which v is constant, there are N(£) 
grid points and on any tj curve ( for which £ is constant ) there 
are N(tj) grid points. Consider any curve £ = £*, which lies 
between the boundary curves £ = £ and £ = £ of the 

min ^ a x 


";rn, 



Fig. 2.2 Eight noded isoparametric element and 
local coordinate system. 



Fig.2.3 Construction of ^-constant lines. 
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domain. Now, on the curve £ = ?\ the grid points are connected by 
a piece-wise continuous curve which is made up of the boundary 
segments of the quadrilateral elements, as shown in Fig. 2.3. 
When Laplace equations are solved, the grid formed by £ = constant 
and 7 ] = constant curves is very smooth i.e., even the slopes of 
the piece-wise segments match at the grid points along the 
coordinate lines. Now we proceed to present the exact procedure 
followed in actual computation. 

The integral statements of the governing PDE'S for grid 
generation can be set up by multiplying eq. (2.6a,b) by test 

. £ tj 

functions V s and V respectively and integrating over the domain 
Q. After integration by parts, the resulting equations for weak 
formulation are: 


{ 


av q a£ 


+ 


av^ a£ 


3x dx 3y 3y 


} 


dx dy 


v? I dl 


0 


(2.7a) 


Q 


3Q 


{ 


dV 0 dv dV 71 djn 

dx dx dy dy 


} 


dx dy 


n 


v 73 JJL- dl 
an 


0 (2.7b) 


3Q 


where n is the coordinate along the outward normal to the boundary 
3Q. According to the Galerkin's approach, the test functions V^, 
V 73 and the solution variables £ and 77 , can be taken to belong to 
the same finite dimensional solution space . For instance, both 
may be represented by piece-wise bilinear polynomial functions. 
Representing the finite dimensional approximations by an over bar, 
eq. (2.7a,b) can be rewritten as: 
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f av^ af av^ af 

| ax ax ay ay 


dx dy 


V? M 


an 


3n 


dl 


0 


( 2 . 8 a) 


f dV V drj + SV 77 37? 

( ax ax ay ay 


dx dy 


V 7 ? 


si 

an 


dl 


0 


( 2 . 8 b) 


n u an 

h h 

Here, the functions V^, V 17 , § and 77 lie in the appropriate finite 
dimensional solution space ( say S ) , while n is an approximation 

h h 

to n, with an h as its boundary. For convenience sake, we denote 
the finite dimensional approximation £ as £ and 77 as 77 
respectively. Let N be the dimension of the finite dimensional 

h 

solution space, S . Expressing £, 77 in terms of basis functions 1 (1 

h 

1 < 1 < # ; we have 

1 h 

N N 

h h 


C = ^ (x,y) 

i =1 


I 


i =1 


#.(x,y) v. 


(2.9) 


As per the standard Galerkin procedure, we assume the test 
functions to be same as the basis functions. Thus, 
vf = (x,y) 

and v? = # (x,y) ,l=si^^ h 

Now, equations (2.8a-b) lead to the following matrix equations 


[ C ] { £ } = {d 1 } 
[ C ] { 77 } = {d 2 } 

where , C = [ C j ] 

is the matrix with entries C = 


(2 . 10a) 
(2 . 10b) 


ay dip. 

dx dx 


f^i a _li } dx dy 
3y 3y > 


1 < i < N , 

h h 

and {d 1 } and {d 2 } are the right hand side vectors whose non-zero 
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entries arise from the Dirichlet or Neumann boundary conditions. 
From eq. (2.8a,b), it is eveident that taking the right hand side 
contribution to be zero ( for £ variable on the r? constant 
boundary and vice versa) , enforces orhtogonality condition at the 
boundary. Since the elements are isoparametric and eight-noded, 
both the physical coordinates (x,y) and the transformed 
coordinates ( £ , 7) ) can be interpolated within each element via 
the expressions 


x 


8 

£>!<*'?> x ± 

i=l 


8 ^ ^ 

y = X ^(X, y) y i (2.11a) 

i = i 




I #,<*■?> 5, 
1 = 1 


8 

7 ) = £ (X ' y) 71 1 ( 2 * llb > 


where (x , y) denote the local coordinates. The relationships 
between lengths can be obtained as follows: 

Setting x = x(x,y) , y a y(x,y) 


dx = dx + ^ dy ; dy = dx + dy 


dx 


(2-12) 


6y 


or 


dx 


' ax 

ax ‘ 


dx 




ax 

ay 




dy 


dy 

ay 


dy 




ax 

ay 





dx dy 

8 a 

i 


dijj 

) — x 
ft Y 


1 = 1 
8 


i = 1 


dX 

Z d\fj ^ 

I 


3>Jj 


X 


1 = 1 
8 


ay 

dip, 


dv 

i=i 1 



dx 


dy 




(2 . 13) 


In a similar fashion , the derivatives can be obtained as 


* fL " 


dX 

ay ' 


■ d_ ■ 


a_ 

3X 

= 

dX 

dx 


ax 

- r j! 

ax 

a 


dX 

dy 


a 

L J 

a 

- ay - 


- ay 

ay - 


- 3y - 


ay - 


( 2 . 14 ) 
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where [J ] is the Jacobian matrix for the transformation which can 
be evaluated using equations (2.13). Making use of (2.14) and 
Gauss-Legendre quadrature formulae , the left hand and right hand 
sides of the matrix equations (2.10) are evaluated for each 
element. The expressions for the basis functions \p are available 
in most of the books on the finite element method (Ciarlet, 1978; 
Oden and Reddy, 1976) . 

Solution procedure : During the actual computations, the frontal 
method is used to solve the matrix system arising from the FE 
approximation of the grid generation equations. It is well known 
that the method adopted for solving the assembled matrix equation 
has a significant bearing on the computer storage requirement and 
execution time when the total number of unknown variables is very 
large. The frontal solution technique is based on the direct 
Gaussian elimination procedure for solving symmetric matrices 
where the leading diagonal is used as pivot. The details of the 
Frontal technique have been discussed ( Irons, 1970; Taylor and 
Hughes, 1981) . 

2.6. Description of the algorithm ; 

In this section, the algorithm for the numerical generation 
of a structured, smooth body-fitting grid inside an arbitrary 
shaped region has been described. For computing the FE solution, 
a guess grid is provided initially using eight-noded, 
isoparametric, quadrilateral elements. The transfinite 
interpolation technique is used for the generation of the initial 
grid. With the initial mesh, the grid generation equations are 
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solved and the values of £ and v are predicted at all the grid 
points of the guess grid. Now, the construction of £ and tj 
constant lines is performed by the use of interpolation in the 
physical domain. The new grid locations x and y are calculated for 
each grid point by determining the intersection points of £ and rj 
constant grid lines. Since the initial guess grid may not possess 
the desired grid properties, the predicted values of £ and 77 (on 
the guess grid ) may not be very accurate. Therefore, the entire 
process of solving the £ -problem and the Tj-problem is repeated 
for one or two iterations with the new grid locations (x,y) 
obtained each time. With each successive iteration, the grid 
manifests more conformity with the given geometry ( except at 
sharp corners ) and the grid lines become smoother. After a few of 
iterations, the grid lines closely approximate ^-constant and 
7 i-constant lines and at this stage, the finally converged grid is 
obtained. These steps of the algorithm are summarized below. 

The corresponding flow chart is given (Fig. 2 . 4 ). 

step. 1 : Define initial data on the given geometry 

Define the boundary of the physical domain fi(x,y). 

Define £ and 77 variation on the boundaries of Q(x,y) and 
Normally, on each 77 - constant boundary, equal increment 
in £ is considered between any two neighbouring grid 
points, and vice versa, 
step. 2 : Calculate initial grid on Q(x,y) 

Obtain the initial grid points x ( i ) y ( 1 ) 
i = with the specified boundary data on 

an, by transfinite interpolation. 
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Step. 3 : Solve for body-fitting co-ordinates ( £ ,v) 

Solve the grid generation equations ( 2 . 6a-b ) 
by finite element method with the boundary conditions 
specified as above in the physical domain. 

Step. 4. Construct ^-constant and ^-constant lines 

Construct £ = £* ( 2 s is N(£)-i) and v = V* (z £ N(t})-i) 

curves by finite element interpolation and 
determine new grid ( x , y ) locations of each 

new new 

,7}j ). Here, ^ and V j are set equal to the values of 
£ and v on the boundary. 

Step. 5. convergence criterion for iterations 

If [ Max |( x ( 1 )- x ( 1 ) | sc (a prescribed 

1 _ . old new , . 

v 1 — i —N 7 tolerance ) 

h 

and 

f Max I ( y old (0- y (0 1 ) - c 

1 sisN ld / 

v h 

then stop the iterative process; 
otherwise go to Step 6 . 

Step. 6 . Perform grid iteration 

Replace the old grid values (x ,y ) by new 

v old i old / x 

grid values (x y ) and go to step 3 . 

new new 

2.7. Important features of the algorithm : 

In this section , the important features of the algorithm 
such as generation of the initial grid, construction of £ and 17 
constant lines etc. have been discussed. 

Algebraic grid generation system 

Generation of the curvilinear coordinate system (£ ,r \ ) can be 
formulated as a problem of finding the cartesian coordinates 
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(x,y) in the interior of the rectangular transformed region from 
their specified values on the boundaries. This can be done by 
interpolation between the specified boundary values of (x,y) , for 
any interior point depending on its location (^ , rj) in the 
transformed domain. The basic problem, therefore, consists of 
finding suitable interpolation functions for generating the 
interior co-ordinates of the physical domain. A large volume of 
literature is available on such algebraic grid generation methods 
and these have been discussed in the surveys by Thompson et al . , 
1982, 1985a; Ghia and Ghia, 1983; Eiseman and Eriebacher, 1987. 

In the present work, transfinite interpolation (Thompson et 
al. , 1985a) has been used to construct the initial grid. The 

bi-directional interpolation is implemented in two unidirectional 
interpolation steps; interpolation is first performed one 
direction (say i i) over the entire domain and then it is applied in 
the other direction (£) , taking due care to match the prescribed 
boundary values of co-ordinates. The procedure is simple to use 
and is inexpensive computationally. It can also be extended easily 
to three-dimensional situations . 

Construction of g and -q -constant lines : 

The construction of £ - and t? -constant lines is done as 

follows. For obtaining the shape of the curve £ = £ , the 

predicted values of the variable £ over the whole domain are 

scanned, to identify the elements on whose boundary £* occurs. If 

for two neighbouring nodes of an element (say j and k) £ < 

* 

- £ , then the £ = £ curve passes between the nodes j and k (Fig. 
2.3). All such elements and the corresponding boundary segments 
are considered for the construction of the £ = £* curve. Now on 
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each of these boundary arcs, it is necessary to find the physical 
co-ordinates (x,y) at the locations where £ = £ occurs. 

Consider a typical elemental boundary arc, on which £* 
occurs. Depending on which side of the element it corresponds to, 
one of the local coordinates ( x or y ) will be constant on the 
arc (see Fig. 2.2). If x is the local coordinate which varies 
along the elemental arc (with y being constant on the arc) , then 
the interpolation specified by equation (2.11a) simplifies to the 
form 

£ = ax 2 +bx +c 

where 

a = ($ t + $ k - 2 £j) / 2 ; b = (€ k - S ,)/2 ; c = Cj 

and i, j, k are the numbers of the nodes which constitute the 
boundary arc. Solving the quadratic equation given by 
a x 2 + b x + c - £* = 0 

the local coordinate x can be determined. Since y has a known 
constant value on the boundary arc, the physical coordinates (x,y) 
of the point at which £ = £*, can be found from equations 
(2.11b) . 

2.8. Results and Discussion : 

The method has been applied for generating grids in simply 
and multiply connected two dimensional regions. Effects of 
boundary spacing and initial grid structure have been studied and 
the findings are discussed below. 

In all the Figs. 2.5-2.11, the intersection of £- and v~ 
constant lines in the final grid are totally different from the 
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initial grid. Within the first one or two iterations, smooth 
grids conforming to the geometry of the problem are obtained, 
irrespective of the initial grid structure. Due to the extrema 
principles of the partial differential equation system in the 
physical plane and the monotone variation of body-fitting 
co-ordinates along the boundary segments, it is clear that 
overlapping grid lines will not occur in the physical domain. 
Further, due to the inherent smoothness of the Laplace operators, 
whatever grid roughness present in the initial grid, disappears 
after the application of the algorithm. 

The grid generated in a simply connected C shaped domain is 
shown in the Fig. 2.5. The domain boundaries are smooth with a 
mild curvature in one direction, while they are non-smooth and 
possess sharp curvatures in the other direction. Starting with a 
rough initial grid, within few iterations, excellent smoothness is 
achieved in the domain interior. Another point to be noted is 
that in the direction of mild curvature , the grid spacing remains 
more or less uniform; in the other direction, grid points move 
away from concave corners towards convex corners. This is a well 
documented property of the elliptic grid generation systems 
(Thompson et al., 1985a). Only a few £ and 7) grid lines have been 
considered for the sake of simplicity, during the generation of 
grids. By incrementing the number of grid lines where required, 
grid smoothness can be improved and grid spacing can also be 
controlled. 

In the Figs. 2.6 and 2.7, grids for some simply connected 
domains have been shown. In one case, a uniform grid spacing on 
the boundary is considered, while in the other, dense 
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Initial grid Final grid 

Fig. 2.5 Numerical grid for the ellipse domain. 
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InitiQl arid Finai grid 

Fig. 2.6 Numerical grid for the simply connected domain : 
Effect of uniform grid spacing. 



Fig. 2.7 Numerical grid for the simply connected domain: 
Effect of dense boundary data . 
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specification of boundary data near the sharp corner of a side has 
been considered. The smoothness of interior grid lines and the 
behaviour near concave and convex boundaries after the application 
of the Laplace solver are evident from the figures. The important 
observation from these figures is that the although the initial 
grids are totally different from each other, the converged grids 
have good correspondence, especially in regions where the boundary 
data remains same . 

The effect of boundary roughness on the grid within an 
elliptical domain is shown in Fig. 2.8. The boundaries of the 
domain have been made rough purposefully, to study the smoothening 
effect of the present algorithm. The transfinite interpolation 
produces a rough initial grid at interior points also. However, 
after the application of the PDE based grid generation approach 
for four iterations, the grid lines become very smooth in the 
interior. As discussed earlier, such grid smoothness is essential 
for obtaining accurate numerical solutions to physical problems. 
It is also seen that the grid spacing is uniform in the interior 
portions, which is a well known property of the Laplace system. 

In Fig. 2.9. the effect of dense boundary grid point data, on the 
O-type grid within an annular region is shown. Due to the convex 
property of the Laplacian operator, the grid lines tend to move 
away from concave corners, towards convex points. The closed t) 
coordinate lines show body-conformity except at sharp corners. The 
£ coordinate lines, however, are distorted due to the ( x-y ) 
symmetry of the problem and the close spacing of the boundary 
points specified at the outer boundary of the domain. 

For uniformly distributed boundary data, both £ and t? lines 
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Numerical grid for the s im P | y connected d 
Effect of non - smooth boundaries 


A B : branch cut 



Initial grid Final grid 


Fig. 2.9 O-type grid for dense boundary data around 
concave corners. 



A B : branch cut 



Initial grid 


Final grid 


Fig. 2.10 O-type grid for uniformly distributed 
boundary data. 


A B : branch cut 



Initial grid Final grid 

Fig. 2.11 0-type grid for arbitrary shaped branch cut 
with a very bad initial guess. 


are smooth and they do not display much distortion, as seen in 
Fig. 2.10. The effects of a very bad initial grid is shown in Fig. 
2.11. The present technique was found to be very stable during the 
grid smoothening process. However, at times, local distortion of 
elements occurred due to the rapid movement of points near concave 
corners. This problem was resolved by under relaxation. It is 
also seen from Fig. 2.11 that a complicated shape of branch cut 
influences the grid only locally- and the effect is forgotten 
further away. 

2.9. Conclusions : 

In this chapter, the Laplace grid generation equations have 
been used to generate body-fitting grids in arbitrary-shaped 
bodies. Finite element technique has been employed to solve the 
grid generation equations in the physical domain itself. The work 
described here has a greate flexibility to account for various 
types of boundary conditions. It is also a rugged procedure 
because of the fact that only linear equations have been used to 
obtain’ the required grid. It must be noted that local refinement 
schemes can be easily implemented by choosing appropriate control 
functions in Poisson grid generation equations and the present 
method can be extended to adaptive grid generation. Extension to 
3-D grid generation is also possible. 
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CHAPTER 3 

FINITE ELEMENT FORMULATION OF FLAME PROPAGATION PROBLEMS AND 
DESCRIPTION OF COMPUTABLE ERROR ESTIMATORS FOR GRID ADAPIION 

3.1. Introduction : 

With the advent of the energy crisis, combustion studies have 
assumed great significance, with a view towards making the 
combustion process fuel- efficient. A majority of the theoretical 
studies on combustion involve numerical simulation, due to the 
inherently non-linear nature of the reaction processes. Central to 
such numerical simulation of combustion phenomena, is the 
capability to track the movement of a flame front. Since the flame 
is a zone of intense heat release and rapid chemical reactions, it 
is characterized by steep gradients of temperature and species 
concentrations. Obviously, adaptive mesh generation is an 
essential feature of any numerical simulation technique for the 
efficient solution of flame propagation problems. 

In the ensuing chapters, algorithms for generating adaptive 
meshes and solving front propagation problems are discussed. In 
the present chapter, general features of the finite element 
formulation and the solution procedure for a one dimensional 
reactive flow system are described. Further, the mathematical 
preliminaries required for defining various measures of computable 
error and the use of these error measures in conjunction with the 
finite element method are presented. 


k 
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3.2. Modelling of flame propagation problem : 

3.2.1. General background : 

The main quantities of interest in any combustion study are 
the global flame properties such as the flame length and the over 
all heat release, the profiles of temperature and species 
concentrations and the flow features like turbulent mixing. 
Detailed theoretical prediction of these quantities is a 
formidable task, due to a variety of reasons. For instance, flames 
typically involve a large number of species and free radicals 
whose rates of production and consumption are governed by several 
simultaneous reactions occurring at different speeds. Another 
major complexity associated with combustion problems is the 
presence of regions in which the gradients are extremely large and 
the transients are very fast. Due to the exponential dependence of 
chemical reaction rate upon temperature, bulk of the combustion 
process is confined to a narrow zone (flame region) , where steep 
gradients of temperature and concentration occur. Tracking down 
the location of the flame region and accurately representing the 
solution gradients in that zone, pose considerable computational 
difficulties. A variety of assumptions and simplifications have 
been invoked in the past during theoretical modelling, to make the 
problem tractable and to obtain solutions of practical importance. 

Traditionally, two distinct approaches have been adopted in 
practice for the calculation of flame structure. Analytical and 
semi-analytical asymptotic methods have been employed to study 
flames described by simple, one-dimensional steady models 
(Buckmaster and Ludford, 1982; Williams, 1985; Kapila and 
Matkowsky, 1980; Ludford and Peters, 1984; Oran et al., 1981). 
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The second approach involves numerical solution of the complex 
system of differential equations for the combustion problem, with 
few simplifications. Today, several powerful numerical methods are 
available for predicting flame structures accurately in 
multi-dimensional problems. Yet, the problems of ignition and 
flame propagation present severe challenges to the available 

numerical methods, because of the small time and space scales 

associated with these processes. 

In recent years, adaptive grid generation techniques have 
been developed for the application of finite difference or finite 
element solution procedures to flame propagation problems 

(Benkhaldoun et al, 1988a, b; Bieterman and Babuska, 1986; Dwyer et 
al., 1980,1982; Dwyer and Sanders, 1983a; Giovangigli and Smooke, 
1989; Larrouturou, 1985, 1986; Ramos, 1985, 1987a-b, 1990a-b; 

Smooke, 1986b; Davis and Flaherty, 1982; Dervieux and Larrouturou, 
1990; Miller and Miller, 1981a). These authors have used the 

concept of distributing the mesh points in a systematic way, so as 
to capture the regions of "large variation" of the solution. Such 
resolution of large gradient regions is important for capturing 
the physical phenomena and also for improving the solution 
accuracy in the flame front zone. The adaptive grid techniques 
reduce computational requirements drastically for a given accuracy 
level and thus make the calculation efficient and feasible. In 
these methods, the location of grid points are determined so as to 
minimize some measure of error in the computed solution. 

Another important feature of the combustion problems is that 
the flame front may be moving with respect to time. Thus, in order 
to capture the flame front, moving mesh methods are needed. For 
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finite difference applications, some authors ( Rai and Anderson, 
1981; Dwyer and Sanders, 1983) have employed the grid speed 
strategy for controlling the motion of grid points. In the finite 
element context, Miller and Miller (1981a, b) and Gelinas et al 
(1981a) , have proposed the movement of grid point locations in 
such a way as to minimize the total residue of the governing 
equations . This moving mesh procedure has been applied for flame 
propagation problems in multi dimensions. 

Benkhaldoun and coauthors (1988a), have proposed two 
different gridding procedures for solving steady laminar flames. 
In their work, the mesh enrichment depends on the reaction rate 
and the adaption is carried out directionally. Also, in a later 
work, Benkhaldoun and Larrouturou (1988a, 1989) have developed a 
dynamic mesh generation procedure based on adaptive refinement and 
derefinement for the propagation of a premixed laminar flame. They 
have used unstructured finite element triangulation and the local 
refinement depends on the second order derivative of reaction rate 
and the diffusive flux. 

In another paper, Benkhaldoun and Larrouturou (1987) have 
combined some features of the finite element and finite difference 
methods for studying the classical thermal diffusion model of 
wrinkled flame front in a gaseous mixture. They have used two 
different adaptive schemes, namely: (i) the 'dynamic rezoning' 
performed at each time step which allows the nodal locations to 
vary smoothly with time, (ii) the 'static rezoning' which is 
performed only at a few time levels during the calculation in 
order to redistribute the mesh points optimally. All the mesh 
points move with the same speed as the flame. The grid velocity is 
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chosen in such 3 way that the int eg ns 1 of each ol the unknowr 
variables on the computational domain roma i nt> constant 01 
equivalently, it is given by the integral of the reaction rate. 

Dervieux et al. (1989b) have reviewed various numerical 
methods for grid adaption in flame propagation problems. Different 
1-D equi-distribution techniques and grid speed methods are 
described. Ramos (1983, 1985, 1990a-b) has used the temperature 
gradient for grid adaption while studying flame propagation by the 
finite element method. In a later work, the author (Ramos, 1991) 
has described adaptive finite element methods based on equi- 
distribution, elliptic grid generation and variational 
formulation. 

In the present work, adaptive grid generation techniques 
which involve a posteriori error measures defined in terms of the 
computed finite element solution are employed. These techniques 
are applied for the study of one-dimensional promixed flame 
propagation and the performance of each technique is assessed by a 
detailed analysis of the error behaviour, solution accuracy and 
the adaptive mesh features. 

3.2.2. Formulation of flame propagation problem : 

The problem considered here is the simplified thermo- 
diffusive model which describes the laminar premixed flame 
propagation in one dimension (Fig. 3.1). A common situation where 
this problem finds application is the combustion of the 
air-gasoline mixture inside a Spark Ignition (SI) engine. In the 
SI engine, ignition is started with the help of the spark plug and 
a flame front then sweeps through the entire combustible mixture 
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A - Unreacted 
mixture zone 

B - Moving flame 
zone 

C - Burnt region 


Fig. 3.1 A SCHEMATIC DIAGRAM OF ONE- DIMENSIONAL FLAME 
PROPAGATION . 
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(see. Fig. 3.2). The propagation of the flame is essentially 
one-dimensional, from one end of the combustion chamber to the 
other. It can be represented graphically in a manner as shown in 
Fig. 3.1. 

The chosen model problem considers the the exponential 
temperature dependence of Arrhenius-type reactions and the 
dif fusional characteristics of more general unsteady combustion 
problems. A channel of unit length, is filled with a cold mixture 
of the reacting species (fuel + oxidizer) in the right 
proportions. A one-step exothermic chemical reaction of the form 

A > B is assumed, where A represents the reactant mixture and B 

represents the combustion products. At time t = 0, the mixture is 
ignited at one end (x = 1.0) by supplying sufficient heat. Each 
calculation is begun by applying a linearly increasing temperature 
at the right-hand boundary (Fig. 3.1). As heat diffuses into the 
mixture, the rate of chemical reaction increases rapidly. The 
result is an acceleration of the reaction zone (flame) which moves 
away from the wall. After the ignition process is completed, a 
self-sustaining flame front develops and propagates into the 
unreacted mixture, until all of material A is converted into B. 
The speed at which the reaction zone moves (flame speed) will 
generally reach a constant value which will be maintained until 
the flame approaches the left side of the channel. This period of 
constant flame speed offers a convenient opportunity for assessing 
the accuracy of various numerical techniques. 

In the present work, the above- described flame propagation 
model has been used for all the adaptive grid studies. Two 
specific cases, one corresponding to a stationary mixtu re and the 
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Intake valve 

Spark plug 

Exhaust valve 

Intake pipe 

Exhaust pipe 

Engine cylinder 

Piston head 

Piston shaft 
(connected to crank ) 

Moving flame front 
(shown as dotted line) 


Fig. 3.2 SCHEMATIC DIAGRAM OF FLAME PROPAGATION INSIDE 
SPARK IGNITION ENGINE . 
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other dealing with flowing mixture , have been considered. The 
model is described by a stiff system of two coupjod nonlinear 
reaction diffusion equations of the form 


li _ ii + v uz 

3 t ax 2 ax 

*-£ _ Ip + v d —P 

3 t ax 2 ax 

with the boundary conditions 


= f(p,T) 

= - f(P,T) 


t > 0, x e (0,1) 

(3.1a) 


3 p( Q/t) _ _ a t ( o, t) 

3 x u “ FIE 


LPd.t) 

3 X 


0 , T (1 , t ) = g (t) 


t > 0 


(3 . lb) 


and the initial conditions 

P(X,0) = 1, T (x, 0) = 0.2, x e [0,1] 
respectively. Here, f(p,T) and g(t) represent the fol .low! ng 
f(p,T) = Q p exp (-K/T) 


g(t) = 


0.2 + t/(2 x 10 4 ) 


t s 2 x 10' 
t a 2 X IQ 


(3.1C) 

with Q = 3.52 x 10 6 and K = 4.0. (31f) 

In these equations, T (= T(x,t)) is the temperature, p ( « p(x,t)) 
is the mass fraction (density) of the reactant mixture, V is the 
flow velocity and the right hand side term f(p,T) is the rate of 
chemical reaction. AH the variables in the problem including the 
spatial and time co-ordinates, have been made dimensionless by 
scaling with the characterstio value of each variable. In the 

Sq ' <3 ' ld) ’ K iS the ^tmensionless activation energy and Q is the 
dimensionless pre exponential factor of the reaction rate. The 


l 
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chosen numerical values of Q and K are typical of a combustion 
reaction. 

These governing equations expressed through a coupled 
non-linear system are stiff due to the existence of a narrow flame 
zone. A few special cases of the general model are considered of 
interest here. If V = 0, the above problem simplifies to that of 
flame propagation in a stationary mixture which can be written as 
follows . 


a t 
a t 

a p 

a t 


a 2 t 
a x 2 

a^p 

a x 2 


= f(P,T) 


= - f(P,T) 


t > 0, x e (0,1) 


(3.2) 


with the same boundary and initial conditions as given by 
(3.1b-f). This problem has been studied in Chapter 4, while the 
general situation of flame propagation in a flowing mixture is 
analysed in Chapter 5. When f(p,T) is replaced by a suitable 
f(x,t), the resulting problem is linear . It can be solved 
analytically and the corresponding solution can be used for 
validating the numerical predictions. In the present study, such a 
linear problem exhibiting a wave front with appropriate boundary 
conditions has been termed as the " Test problem " while, the 
general non-linear reaction diffusion equations ( 3.1 a-f ) are 
referred to as the " Flame problem " . The test problem case has 
been used to validate the numerical solutions of problems 
considered in Chapter 4 and Chapter 5. 


3.2.3. Salient features of the flame propagation problem : 

In the reactive diffusive model presented above, the chemical 
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reaction rate exhibits a strong non-linear dependence with respect 

to temperature. In fact, according to the " Arrhenius kinetics", 

the rate of reaction depends exponentially upon the mixture 

temperature. One of the main difficulties in simulating llame 

propagation described by such reaction kinetics is the existence 

of disparate characteristic length scales. The stiffness of the 

problem appears essentially in the form of a narrow zone of sharp 

gradients around the flame front, dividing two regions in which 

the variables are almost constant. Thus, accurate representation 

of the inner structure of the thin flame plays a vital role in the 

overall simulation of the problem. Other aspects of interest in 

the study of the flame propagation problem are the predictions of 

the flame speed and the ignition transients. The critical step in 

the entire solution procedure is the accurate resolution of the 

flame front region at each time instant. 

• 

3.2.4. Computational requirements of flame propagation studies : 

One of the major difficulties encountered _ during the 
simulation of moving front problems is the lack of information 
concerning the gradients of the dependent variable in different 
parts of the solution domain. Without this vital information, high 
gradient regions can not be resolved in a computationally 
efficient manner. Also, if a uniform grid is employed, a very 
large number of grid points will be needed to produce accurate 
results. However, truncation error on such a non-optimal grid will 
corrupt the computational results, because of the unnecessary 
usage of a small grid size in most of the solution domain. 
Besides, the requirement of obtaining a solution of given accuracy 
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with the minimum total cost and. computational effort can not be 
met with the uniform grid. Thus, the chosen numerical method must 
be in a position to bunch grid points in high gradient regions and 
spread them out in regions where the solution varies weakly. It 
is expected to reduce the computational effort significantly 
compared to that for the uniform grid, thereby making the 
calculations efficient and feasible. These needs call for the use 
of a powerful adaptive grid procedure. 

3.3. Mathematical formulation and finite element solution 

procedure : 

The need for adaptive mesh generation during the solution of 
the front propagation problems, was amply elaborated in the 
previous sections. Adaptive grid refinement can be applied to the 
front propagation problems by employing various numerical 
techniques such as finite difference or finite element methods. 
The immediate aim of our work is application of adaptive finite 
element solution procedure to the problems as described earlier in 
which, the adaptive grid refinement is based on some computable 
error-, measured in a properly chosen norm associated with the 
problem. This, in turn, requires the definition of a proper a 
posteriori measure of error so that the grid refinement is 
carried out at locations of large error, either in the solution or 
in the gradients of solution. For characterizing the predicted 
solution at each stage of grid adaption and for identifying the 
error measures, functional analysis tools have been used. In the 
ensuing sections, some important mathematical concepts based on 
functional analysis are presented which are useful for the grid 
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adaption as well as the finite element solution of flame 
propagation problems. 


3.3.1. Mathematical Preliminaries : 

In this section, we introduce some notations and definitions 

which will be used throughout this work. All functions treated 

, Q J f 

here are real-valued functions. For each integer j^O, D f = — - . 

Ox J 

Let Q be a domain of variation of x independent of time, such that 

0 < x < 1. On the space of square integrable functions on Q, the 
2 

L -inner product and norm are denoted by 

<u,v> = J uv dx 
0 


and, 


II u II 2 = < u,u > 

respectively. The support of a function f is defined as 

Supp f = { x : f (x) £ 0 } . 

Let Q be the closure of Q; 8 Q be the boundary of Q and 

Sfi = 30 u df2 + . Also, 0 = n U SO. For any non-negative integer a, 

(X 

C (Q) denotes the space of functions with continuous derivatives 
up to and including order a in Q. As usual, C™(Q) denotes all 
C m (Q) functions having compact support in Q, while c"(Q) is the 

space of all infinitely differentiable functions with compact 
support in Q. 

For each non-negative integer m, the Sobolov space of order m 

on n denoted by H*(Q) is defined as the closure of C m (Q) with 

respect to the norm 


II u II 


H m cn) 



D J u (x, t) | 2 dx. 
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Equivalently, H (Q) is the space of all functions whose 


distributional derivatives of order less than or equal to m, are 
in L (Q) . Clearly L (Q) = H°(n) . H m (Q) is a Hilbert space with 


the inner product defined by 


<U ' V> H-(Q) 


• r 

* E J 

*T\ « 


D J u D J v dx 


For convenience sake, we use the notation 

* -£. * - & *, = I* 

x OX xx 2 t d t 

oX 

We shall also employ the convention that the variables 
i,j,l,m,n,j and 1 take integer values. The the total number of 
partial differential equations under consideration is denoted by 
npde . On 0, we consider the family of partitions 9 


I = A(t) = j 0 = X q <X J < 


< X * < X * = 1 € 9 

l-i l f 


where A(t) is a space grid of n. The following symbols are used to 
denote an element (0^), the boundary of space domain (50), step 
length (h^), maximum step length (h) and minimum step length (h ). 


0 = Q (A) 

j j ; 


X€[R:x <x<x 

j j+i 


j = o,i, . . . , 1 -1 


so = an" u an = { o,i } 


h = h (A) 


\ h (A) 


Oij^l-l j 


h = h (A) = x - x ; 
j j ‘ j+i j 

h = h (A) = J 1 " * h (A) 

- v ' oijil-l j 


9 is said to be a K(X) -regular family if constants k i and X > 0 
such that h (A) a: X h*(A) for j = o, i , . . . , i *- i . For any vector 
y e R 1 with components y , y j7 y g , . . . .y l , define norm by 

IMI ro = Max .K) 

0^1 < 1 
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3.3.2. Basic formulation in a general set up : 

The reaction diffusion equations (3.2, 3.1b -3.lt) for the 
flame propagation in a stationary mixture can be incorporated into 
a general system of non-linear parabolic problems with the 
appropriate boundary conditions. In this section, we describe the 
finite element solution procedure for these equations which 
involves weak formulation, Galerkin procedure and time 
discretization scheme. 


(a). Description of the general problem : 

The reaction diffusion equations (3.2, and 3. lb-3. If) for 
the flame propagation in a stationary mixture can be written in 
the more general form : 


(u,) t - 




f t (x,t,u) 


(x,t) e Q X ( 0 , t ] , 

1 = 1 , . . . , N P I! K 

(3.3a) 


with the boundary conditions 

V x) + V*) (Vx = 9 x (x,t:) 


(x, t) e c?fi x (0, t ] , 

1 = 1 , . . . , NPDE 

(3.3b) 


and initial conditions 

Uj (x, 0) = u^x), X € Q (x,t) € Q X (0,t f 3, 

1 = 1 , , NPDE 

(3 . 3C) 

Here u (= u(x,t)) is a vector which consists of entries 

u i (x 't), u 2 (x,t) , u NpDE (x-,t) . 

For the flame propagation studies in the present work, the value 
of npde is 2. Further, the solutions u^x.t) and u p (x,t) 
represent the mixture temperature and density respectively. 


In 
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the following discussion, we may refer to the solution vector by 
its component u^ or simply by u. Also, it is to be noted that u 
appears as a vector in the non-linear reaction term f (x,t,u) of 
the above equation. The coefficients a i (x) and /3 } (x) determine 
the type of boundary conditions and S <x) is the initial condition 
at time t = 0 for each u x . For x € £2, we consider a (x) ^ 0, with 
a^x) > a Q > 0 for atleast one index u Further, it is assumed that 
(3 i (S£2 ) ^ 0; /3 j (d£2 + ) £ 0 and a^x) £ 0 for x e 3£2, i = i, npde. 
We also assume that with reasonable problem data, the solution 
u t (x,t) of the system of equations (3.2, 3. lb-3. If) exists and a 
smooth mapping of [0,1] into the Hilbert space 


H* pde (£2) = 


{v} : v,(v) e L (£2) ; i = i, npde 

1 I ’ 1=1, NPDE 1 v r x 2 ' 1 ' 


Let us define the seminorm 


rl f UL 1 / 

{ X < M v ,>» <v*> } 


v e H 


which is a norm on h” pde (£2) provided appropriate conditions on 

a a and g ( i s i s npde ) hold. 

1,1, i 


(b). Weak formulation : 

Now the weak solution u i (x,t) (= u^ of (3.3a-c) is defined 

as follows : Find a function 

u : [0,t f ] > H x (£2) i = !/•••/ npde . 

such that 

o^Uj^t) + /3 i (x)(u i ) x = g^Xjt) for i = i,...,npde on x 6 
is satisfied and u i (x,t) satisfies the equations 
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<lt a ' < -' t) - v + < a , s i|<x " t) 4 Vi > 

+ a x (x) a % (x) Gj (x,t)v i 


/3, (*) 


an 


an 


+ 


a (x) i an 

< f x (X,t,U),V ) > + g i (X,t)V i 


+ 


an 


t e (0, t ] , 1 = 1 , - . . , NPDE . 


(3.4a) 


and < Sj (x, 0) , v i > = < u(x), >, t e (0,t f ], i = i , . . . , npuk. 

(3.4b) 

for all 1^(0). In fact (3.3 a-c) to (3.4 a-b) are equivalent 

if u i (x,t) has sufficient regularity. Let us now turn to the 
formulation of Galerkin method in space variable called the 
semidiscrete approximation for the solution 

u(x,t) = ^(Xjt), u 2 (x,t), u NpDp (x,t) 

of the equations (3.4 a-c). 


(c). Semi discrete approximation : 

We describe the construction of the finite dimensional 
solution space in which the approximation for the solution u i (x,t) 
i n H p (1 < p < <» ) lies. Generally the finite dimensional solution 
space, considered to approximate u (x,t) in the Galerkin 
procedure, belongs to a regular family and satisfies certain 
inverse properties. In the finite element literature, such a 
class of finite dimensional sub spaces which depend on the mesh 
size (h) is referred to as S^' ^ family. We recall here the 
definitions of S h ' families as described in Babuska (1972) and 
Oden and Reddy (1976) . 
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Construction of finite dimensional solution space : 

Let ri 1 be a given integer and k = ( k 0 > k i k* ) be a 

vector of integers corresponding to the 1 sub— intervals such 
that 

2r, j = 0 

Then S^ ,k ( = S^’ k (A) ) is defined as the space of (r-1) times 

continuously differentiable functions on Q for which the 

restriction to (j = 0 , 1 , 2 , — ,1 - 1 ) is a polynomial of degree 

k -1 (Ciarlet, 1978) . In fact, for any u e S r,k (A) , let (u ) 

J 1 h 1 j 

be the restriction of u j to j = 0 , 1 , 2 , , i*-i. Then for 

kj= 2r, the polynomial ( i^)^ is defined by its value and those of 

its (r-1) first derivatives at the end points of Q . For k i 2r, 

j j 

( Uj)^ is a linear combination of polynomials for which the value 
and those of the (r-1) first derivatives are prescribed at the end 
points of Q^and k^ polynomials which together with their (r-1) 
first derivatives vanish at both these end points. Clearly, we 
have S r,k (A) c H r (Q) for any 1 < p < « . Then a continuous time 

h p 

Galerkin approximation of the solution 

U t (x,t) = fu x (x,t), U 2 (x,t), U NpDE (x,t)) (3.5) 

is defined by a function 
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with suitably chosen initial conditions 

< U (x,o) , v > = < Gj (x) , V l >, t € (0,t f ] , and 1 = i , . . . ,npde. 

(3.6b) 

In the above expressions, 6 i (x) e S^ ,k (A) is the L ? projection of 
U (x) in S r,k (A) and f (x,t,U(x,t)) is an interpolation of 
f i (x,t,U(x,t)) . 

Expression in terms of the basis : 

Let \p (x) <= S r,k for each j = 1, 2,...,N, be a basis for 

1 t J h u 

r k 

the finite dimensional solution space S ’ for each i = i, npde 

* h 

. r k ... 

where dim (S’ ) = N . Then the Galerkm approximation U (x,t) 

h h l 

can be written as 


N 

h 


u 

1 

(x,t) = l U [t] xfj 

J = i 

. (x) , 

» J 

t 6 

(0,t f ] , 1 = 1, . . . , NPDE . 







N 

(3.7) 

Here, 

the 

vector Ujft] 

- K 

> [tJ ' 

I h 

V represents the solution 

' i=i 

vector 

for 

each l = i , . . . 

. . , NPDE 

and 

in the flame problem 

this 

vector 

is 

equivalent to 

temperature 

and density at each 

nodal 


value of the mesh A. Substituting the above expression of U^Xjt) 
and taking 

= V i for 1 = i , . . . , NPDE 

in the equation ( 3 . 6a-b ) , we obtain an initial value problem for 
the system of N h ordinary differential equations for each partial 
differential equation. Define 

M , = [( M j),,] as the matrix with entries (H ) = <\b , xb > 

i i i j v r ij v i i ' M j 

K i = J as the matrix with entries (K ) = <a \7i// ,Vi h > 

1 1 1 J i ij i v n ' v i j 
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F i = * F ij } is the ri< 3ht-hand side vector with entries 
F = <f, (x, t,U[t] ) , \Jj > 

r i = as the vector w ith entries = U [0], which are the 

components of initial approximation ft (x,0) and the dimension of 
all these quantities equal to N , the dimension of S r,k (A). Also, 

h h 

for all these variables, the subscript ' i ' varies from i to npde. 
Employing the above matrix notation, equation (3.6a-b) can be 
written in an equivalent form as : 

[mJ{ | + [ k] {u[t]| = (U[t] ) for 1 = 1 ,..., NPDE 

(3.8a) 

with the initial conditions 

| 1 for l = l . . . , npde (3.8b) 

The matrix J is positive definite V t ^ 0 and the ODE system 

has a unique local solution for t ^ 0 and i = i,...,npde. 

3.3.3. Method of solution : 

The system of ordinary differential equations (3.8a-b) in 
' U [ t ] j = 1 , . . . N 1 ; i = i,..., npde 

1 , j h 

are nonlinear as F (U[t]) is nonlinear in U[t]. Thus, the 
determination of U i (x,t) amounts to solving the system of 

differential equations in U [t] with appropriate initial 

1 > J 

conditions. These equations can be solved by a variety of 

numerical methods available for the solution of a system of ODE, 
one of which is to discretize in time and generate a set of 

simultaneous algebraic equations. Several researchers have used 
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the integration packages which are designed to solve specific 
classes of ordinary differential equations accurately and 
efficiently (see Hindmarsh and Gear, 1979; Axelsson, 1976; 
Bieterman and Babuska, 1986; Benkhaldoun and Larrouturou, 1987). 
Many ideas of these solvers have been utilized in the context of 
evolution problems for PDE's and particularly in conjunction with 
finite element methods. As the resulting difference equations are 
nonlinear, some authors (Davis and Flaherty, 1982; Ramos, 1985) 
have adopted the incremental Newton-Raphson method. It is well 
known that the discrete representation of the ODE operators have 
certain properties that lead to restrictions on time step and the 
element length, to reduce phase errors which generate 
oscillations. The choice of a particular scheme for solving these 
ODE systems, therefore, depends upon the step size restrictions 
and the accuracy of the method. 

In the present study, the Crank-Nicolson-Galerki n method has 
been employed to discretize the problem with respect to the time 
variable. This procedure is second order accurate in terms of the 
time step At. The resulting non-linear algebraic equations are 
solved by an iterative method. The iterative method converges 
free of oscillations if adequately small time steps At ( of order 
10 5 to 10 6 ) are chosen. 

(a). Time discretization : 

Let At be the time step and U*(,,t) be the approximate 
solution m s^’ k (A 1 ) of DJ^t) at t = t ( = i At (i s i s npde) 
Here, the finite dimensional solution space S*' ,k (A l ) is 
constructed at the time t t of the mesh A with dim (S^’ k (A 1 )) = 
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The semi— discrete equation is then discretized with respect to 
the time variable in a symmetric fashion around the point 
t i _ i/2 = ( 1 - 1 / 2 ) At producing a second order accurate method in 

time. More precisely, we define u|(.,t) in S r,k (A 1 ) recursively by 

1 h 


i-1 / 2 

a u, (.,t) 


<K l ' ' *i> + < £ 


1 Sx 


ujf-t) + u;- x ( .,t) 

' ax l > 


1 < fj +~f;- 1 (. / t,U i - 1 (.,t)), 

for ip i in S r ’ k ( A 1 ) ; t e [t ,t ], 1 = i,...,npde; 


(3.9) 


with appropriately chosen boundary conditions. The expression 
for the term 1/2 (.,t) can be derived as follows: 


Uj(.,t) = U|" 1/2 (.,t) + ~ Uj" 1/2 (.,t) + 0 ( (At) 2 ) 

Uj _1 (.,t) = u;- 1/2 (.,t) - ~ U X_1/2 (.,t) + 0 ( (At) 2 ) 

U|" 1/2 (.,t) = Uj _1 (.,t) + — |t U;- X (.,t) + 0 ( (At) 2 ) 

u;- W 2 (.,t> = uj(.,tj - It u ! ( -' t} + 0( (At)2 } 


(3.10a) 


(3.10b) 


(3 . 10c) 


( 3 . lOd) 


From these equations, one can write 

1 1-1 


1 - 1/2 


au (.,t) _ a 

nl. ■*> 


at 


at 


u x (. ,t) + U(. ,t) 


At 


Uj(.,t) 


uj _1 (. ,t) 


(3.11) 


Thus the final discretized equations can be written as 


< uj(. ,t) - U x_1 ( . , t) , + At < /a i ~ 


a f u|(.,t) + u;- x (.,t) 


N a ip 
' 1E 1 > 
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= At / f ) (• ,t,U l (. ,t) ) + f | _t (• ,t,U l " 1 ( - , t.) ) , i // t ) , 
for 0 in S r,k (A 1 ) ; i = i , • . - , urn:? t * [ t s _ x j 

(3.12) 

It is noted that in the above system of equations (3.11) the 
boundary conditions have been incorporated in the choice of iji 's. 
Thus the matrix equations for the time discretized form (3.12) can 
now be written as follows : 


' M, + K, ] { u; ( t] } 


At 

M - — K 

l 2 




= At 
2 


|fJ (U 1 [t])| + At (U 1_1 [t]) 


t e [ t , t ] ; i ~ i , . 
L i - i ' i J ' 


N 


where u| [t] = u| ^ } [t]^ (x) ; t e [t^^tj; i 


, Ni'ilK 

(3.13) 

, NI’HK 


J=1 


Since the solution vector, namely, 

1 [t]| for each i = i , * * * , mm:; 

is known, the second term on both sides of the equation ( 3 . 13 ). is 
known. Thus the matrix system can be rewritten as given below for 
each t e [ t^, tj. 

[ A l] { u ! Ct]} = { G I < u' [t] >} + | hJ-‘ (U 1 - 1 [t]| 

i = npde; t e ft ,t 1 

L s -1 ' 1 1 


(3.14) 
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where 


[ A l] ■ [ M , + “ K , ] 

{g[ (U 1 [t] ) )- = it|pj (U 1 [t])| 

| H ;- I (u , - I [t])} - [ M , - k,] {u;-*[t])} + At |F;-‘(u i - i [ t:)} 

(b). Iterative method : 

The above equations (3.14) can be solved by the following 
iterative method. The main aim is to obtain the finite element 
solution at t = t by solving the matrix system of equations 

(3.14) with the help of the finite element solution available at 
t = t ; _ i .Since the unknown variable appears in both sides of the 
equations (3.14), the corresponding equations are nonlinear. A 

proper initial guess solution is very important in order to obtain 
the convergence of iterative method. For the starting iteration, 
the initial guess solution for the evaluation of source vector at 
t = t is provided by the finite element solution available or 
computed at previous t = i . But, at the subsequent iterations, 
the freshly computed solution at t = t is used in the source 
vector at t = t . Now, the right hand side of the matrix 

equations (3.14) is evaluated and then the matrix system (3.14) is 
solved by direct methods. Finally, the iterative procedure is 

repeated and stopped after the maximum difference of successive 
solution vectors becomes less than a given tolerance ( e ) which 
is of the order (10~ 6 ) . The convergence of the iterative method 
poses no problems since small time steps have been considered. 
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This iterative procedure is explained by the following steps of 

A 

the algorithm. In the algorithm, we denote U[t ) as the staring 
guess solution for the iterative method, U[tJ the computed 
solution, M the total number of iterations performed and 


A i 

u [t] 


M 


as the solution at M th iteration respectively. 


Step :1 Choose the initial guess solution at t = t (1 


1 £ NI’DK). 


i . e . , choose 


r a i i M 
( V^J - 


uj [t] = uj ' 1 [t] (1 


Nt’DK ) 


Step 

Step 


where M is the iteration count set equal to 0 initially. 
:2 Increase the iterative counter M by 1. 

:3 Calculate the right hand side vectors at t = t 


i.e., {g[ (u‘ 


[t] 


M 


l at t=t t and j h| _1 (U 1 ” 1 [ t J ) ) at t=t 


i- 1 


step :4 Solve the matrix system (3.14) for each i (i £ i nim>k) 


Step 


and denote the solution by u| [t] 

5 Convergence check 

M , , . \ M+l i 


ujtj 


MH 

at t - t 


i 


If 


r a i i r a i i 
u t [t] - U [t] 

V /V > 


e for each wnm: ) 


M 


Step 


f A « \ 

then accept jl^ft) = U‘ [t] as the computed solution at 
t = t t and stop the iterative method. 

:6 Perform iterations 

If the convergence criterion (step 5) fails then 

, f A t , J+l 

replace ( U [t]j byj U l [t ]j and go to step 2. 


3.4. Description of error measures used for grid adaption 
in the finite element method : 

In this section, we describe the various error measures used 
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for the grid adaption in the literature. 

3.4.1. History of a posteriori error estimates used for grid 
adaption : 

Several researchers (Babuska and Rheinboldt, 1978a, 
1979b, 1979c, 1981a, 1982a; Bank, 1983, 1986; Bieterman and 
Babuska, 1982a, 1982b; Carey and Humphery, 1981; Eriksson and 
Johnson 1987, 1988; Gago et al., 1983; Oden et al., 1986a, 1986c; 
Oden and Demkowicz, 1986b; Szymczak and Babuska, 1984; Zhu and 
Zienkiewicz, 1991) have used the a posteriori error estimates 
based on residuals as a basis for the local enrichment of the 
solution. The analysis of error estimates and refinement 
strategies have been greatly aided by the the recent developments 
in the mathematical theory of finite element method. For this 
purpose the concept of error estimators and indicators have been 
developed, the estimators giving a relatively accurate measure of 
the error in a given mesh and the indicators a measure of change 
in the energy associated with the introduction of new variables. 
These a posteriori error estimates provide a basis for the design 
of adaptive finite element solvers. 

For finite element grid adaptions, the h and h-p refinement 
techniques have been widely used by several authors ( Babuska and 
Rheinboldt, 1982a; Babuska et al., 1984; Demkowicz et al., 1989; 
Oden and Demkowicz, 1986c; Oden et al., 1989; Szabo, 1986; Gago et 
al., 1983, Zienkiewicz et al., 1989) in which the convergence 
properties which can be improved, there by reducing the error in 
the finite element solution. The practical implementation of the 
optimal h-p refinement strategy is of difficult and has not yet 
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dense meshing and high p degree approximation procedure;-, which are 
not simply adapted to practice. Most of the general ideas on the 
error based refinement strategies for finite element app] ieations 
have been presented in the survey paper by Babuska and Rho j. nbo Ldt , 
1982a; Oden and Demkowicz, 1986c. 

Babuska and Rheinboldt (1978a, 1981a) have derived computable 
a posteriori error estimates for finite element solutions in an 
asymptotic form as h — » 0 where, h measures the size of the 
elements and the grid refinement depends upon the a posteriori 
error estimate. Further, the same authors have also derived the 
error estimates in order to achieve an optimal mesh in which all 
the error indicators are nearly equal. 

These ideas have been extended to parabolic problems 
involving a linear self adjoint operator of the second order by 
several authors ( Bieterman and Babuska, 1982a, b; Kriksson and 
Johnson, 1988; Lohner, 1987; Oden et al., 1986a). Bieterman and 
Babuska (1982a, b) have derived computable a posteriori error 
estimates of the space discretization error in the finite element 
method. The effectiveness of the error estimation is related to 
the regularity of the solution, mesh family type and asymptotic 
range for the mesh size. The theory presented in their work, 
provides the basis for the analysis and adaptive construction of 
time dependent meshes in spatial region. With the help of such 
spatial error and discretization error estimates, the grid spacing 
can be selected so as to equidistribute the local error and, grid 

refinement can be continued until the total error meets a certain 
tolerance criterion. 

The availability of precisely defined reliable 


error 
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estimators appears to be central to the design of adaptive 
procedure, but for the nonlinear problems like flame problem, the 
error estimates have not yet been developed although such a 
estimates are important from the computational point of view. 

3.4.2. Description of a posteriori error estimates 

In this section, we describe the important error measures 
for time dependent problems which are useful in the present work, 
(a). Description of the model problem : 

Following Bieterman and Babuska (1982a, 1982b) we describe 

the following model problem and the various error measures used 
for the linear parabolic problems. Consider the following model 
problem : 

|^u(x,t) ( a ( X )§5T~ ) + b(x)u(x,t) = f (x, t) 

(x , t ) e I x (0,t f ) (3.15a) 

with the boundary conditions 

u(0, t) = 0 = u ( 1 , t) ; t € (0,t ) (3.15b) 

and initial conditions 

u(0,x) = u q (x) ; x e I (3.15c) 

Here, I = (0,1) is an unit open interval in 1R 1 . Further, let a(x) 
and b(x) be sufficiently smooth functions on I for which 

0 s j| < a (x) ^ a < oo ; 0 ^ b(x) ^ b < oo VxeT 

Let T be a family of mesh partitions of I and 

A = i 0 = x <x< < x„ .,= li € ? 

( o 1 n ( A) J 

where N(A) is the dimension of corresponding finite dimensional 
solution space in which the solution lies. Let h= x - x and 

= [X-s/ x ]. With reasonably smooth data for u and f, we 

j L 3 ' j + i J J o . 
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assume that the solution of equations ( 3 . 15a~c ) satisfying 

certain regularity properties. The details of the regularity class 
can be seen in the references (Babuska and Rheinboldt, 1978; 
Bieterman and Babuska, 1982a ,b) . 

Now for each A e T, the finite element approximation is 
obtained in a routine way as explained in the section ( 3.3.3 ). 
Let u(x,t) be the exact solution and U(x,t) be the computed 
solution of the system of equations (3.15a-c) respectively. 

The a posteriori error estimates of the elliptic problem play 
a crucial role for the development of the necessary theory for the 
parabolic problems. Here, the model problem for the elliptic case 
(see Babuska and Rheinboldt, 1978a) is deduced by dropping the 
time derivative term in the equation (3.15a ) given by 

“ §~ x ( a ( x )§~^ u ( x ) ] + b (x) u (x) = f ( x ) ; x e I 

(3.16) 

with the appropriate boundary conditions. 

(b). h refinement technique : 

A posteriori error estimate for elliptic problems : 

The following asymptotic a posteriori error estimate for 
finite element solutions has been derived by Babuska and 
Rheinboldt (1978a) for a class of linear elliptic equations of 
type (3.16) . 

N(A)-1 

g/z 

ll e ( x ) II E * c £ [l+0(h)J as h -•* 0 

i J=0 

( 3 . 17 ) 

where, e(x) = u(x) - U(x), u, U being the exact and computed 
solution of equation (3.16) and the energy norm of the exact 
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error (||e|| ) is given by 


1 e ( x ) || g =^a(x)|e(x), de(x)y + ^b(x)e(x) , e(x)^ 

(3 . 18) 

Here, c is a constant which depends upon the degree of 
approximation in the finite element method. In the equation 
(3.16), the computable quantity 7) is called the error indicator 
and is expressed by : 


,2 x 
h ,j + 1 


3. r v 

TT 2 X i 


dx , j = 0 , 1 , 2 ,... ,n(A)-i 


(3.19) 


where the elemental residual for the element Q is given by 


r = r (x) 
j j v ' 


d_ 
d X 




+ b(x)U(x) - f (x) , 


X € [X , X ] , 
L j ' j+1 J ' 


- ( X + X ) 

2 v j j+l' 


a = a | 

and N(A) is the number of nodes in the domain. The a posteriori 
error estimate in the equation (3.17), namely 

' N(A> '^ ) 1/2 

Z V J (3.20) 

J=o 


provides an easily computable bound for the total solution error 
under the energy norm. Further, the error indicator tj provides a 
measure of the error based on the local elemental residue r^ of 
the j th element of the mesh. The coefficient 1/n 2 which can be 
replace by 1/12 in the equation (3.19 ) has been obtained by an 
eigen value analysis of the associated operator satisfied by 
finite element residue of the approximate solution. 

A posteriori error estimate for parabolic problems : 

Using this error analysis for the elliptic problems. 
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Bieterman and Babuska (1982a, b) proposed an a posteriori error 
estimate for the time dependent spatial meshes and presented the 
theory for the estimation of energy norm of the space 
discretization error at each time. In order to assess the 
performance of the a posteriori error in actual practice the 
authors have used various notions such as relative error, true 
error and effectivity ratio. We define some of the error measures 
which are used in conjunction with the time dependent problems. 

The energy norm of the exact error is defined by 

||e(.,t)||* = / a (x) 8 e(x,t), 5 e(x,t)\ + / b(x)e(x,t), c(x,t) \ 

h ' ax ax ' ' ' 

The gradient norm of the exact error is defined by 
III e (-,t) HI = | | a (x) § x (e(x,t) ) 2 dx j 



r j(x,t) = |_u(x,t) - ~ a(x)|^ (x) J + b (x)U (X, t) - f(x,t) 

(x,t) e QjX (0,t f ) . 

It is computable for each t > 0 on each with the help of the 
computed solution U(x,t) of the equations ( 3 . 15a-c ). The 
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Relative error E REL (t) is defined as the ratio of the computed 
gradient to the exact gradient in some properly chosen norm 
associated with the problem and is given by 

E REL (t) = l|e(.,t)!/|U(.,t)l (3.22) 

Effectivity ratio 9 ( t) is given by 

0(t) = 8 (t) / ||Je ( . , t)||| (3.23) 

and it indicates relative growth of a posteriori error and 
gradient norms of the error. Bieterman and Babuska (1982a, b) have 
shown the existence of asymptotic estimates with respect to 
spatial mesh size for the energy norm of the error ( fle ( . , t) || ) . 
The mathematical theory presented in their work provides the basis 
for the estimation of ||e(.,t)|| of the Galerkin approximation of 
time dependent problems. Further, they have carried out an 
asymptotic analysis for the estimation of energy norm and have 
shown that for higher terms in the space mesh size 

||e(.,t)|| E _ HI e ( - , t ) If 

The latter error estimate ( |||e(.,t)j|| ) of the above equation is 

related to the computed solution gradients. Since the exact 
solution is not available, the a posteriori error estimator 8( t) 
has been used in place of ||| e(.,t) ||| for the adaptive procedures. 
In order to assess the performance of 8 ( . ) , which depends on local 
elemental error indicators, the concept of effectivity ratio 0(t) 
has been developed by these authors. Further, they have studied in 
detail the behaviour of 8(.) and the realtive error E rel (.). From 
the above definition of 8( t) , it can be said that, g(t) is an 
effective estimator if there exists a reasonable constant 'C 
(depending on the admissible solution class for Eqs . (3 . 15a-c) and 
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the admissible class of mesh partitions but not on the magnitude 

of the data U Q> f, and g) such that 

l/0(t) as c v t > 0 , 

In practice, it is desirable that 0(t) be close to unity in 
order that the computed error is a proper measure of the true 
error. Therefore, as the partitions of (0,1) are refined 

0 (t) » 1. The quality of the estimator f.'(t) was theoretically 

analyzed for the uncoupled partial differential equations. It is 
to be noted that relative error expression signifies how well the 
gradients are approximated. This is important for front 
propagation problems due to the presence of high gradient regions. 
The authors have also applied the h refinement technique for the 
reaction diffusion type of problems. 

(c). h-p refinement: 

Extending the work of Babuska and Rheinboldt. (1978a), Kelly 
et al (1983) and Gago et al. (1983) have suggested an error 
estimator for higher order elements for linear elliptic equations 
°f type (3.15). They have used the super-convergence results at 
nodes in their assumptions of exact computed values. (i.e. 
e(x.) = 0 for each grid point x i e I). Here, the local error 
expression for the a posteriori error estimate is given by 




12 a p 2 


r J+1 2 . . 

i h (x) 


dx 


(3.24) 


where p^ is the local degree of approximation used in the element 

fij and r^(x) is the local elemental residue of the governing 
equations in . 
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3.4.3. A posteriori error estimates used in the present work : 

In the present work, we have used the error estimates 
proposed by the Babuska and co-researchers (Bieterman and Babuska, 
1982a, b; Kelly et al., 1983; Gago et al., 1983) for the 
computational simulation of both the Test problem and the Flame 
propagation problem in a stationary mixture, by the finite element 
method. Further, the error estimates have been used for the h and 
h-p refinement techniques in Chapter 4. Since the flame 
propagation problem is non-linear and time dependent, no 
mathematical theory has been been developed for the validation of 
numerical simulation of the present work. However, the a 
posteriori error estimates for both refinement methods have been 
tested using the analytical solution of the test problem. In 
addition to the error estimates of Babuska and coworkers, a novel 
self correcting error estimator strategy has been proposed and 
used in the simulation of flame propagation in a flowing mixture 
in Chapter 5. The grid generation algorithms and the predicted 
results are presented in Chapters 4 and 5. 
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CHAPTER 4 

ADAPTIVE GRID TECHNIQUES FOR FLAME PROPAGATION PROBLEM 

IN A STATIONARY MIXTURE. 


4.1. Introduction : 

This Chapter' describes the h and h-p version adaptive mesh 
strategies and their application to the flame propagation and the 
test problems. A posteriori error criteria based on residual 
error have been used and tested numerically for the time dependent 
meshes. A statistical approach for the mesh refinement in h and 
h-p version has been explained in detail. The decisions for 
applying h-p refinement strategy and the conversion of higher 
degree to linear degree approximation in the h-p version have been 
explained. The effectiveness of the a posteriori error estimator 
has been studied by comparing the true and the computed errors at 
various times for the test problem in both the refinement methods. 
The conclusions drawn from the behaviour of various types of 
errors have been discussed. An attempt has been made to construct 
computationally inexpensive meshes for uniform front motion by the 
grid velocity concept. The performance of the h, h-p and grid 
speed schemes have been compared for both the problems. 

4.2. History of error based adaption strategies : 

The adaptive construction of meshes for the numerical 
solution of time-dependent partial differential equations has 
attracted wide interest in the recent past. Particularly appealing 
are the methods which involve the equi-distribution of some 
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measure of error, since they provide a well-graded mesh along with 
control over the solution error. Several schemes have been 
proposed based on the finite difference and the finite element 
approaches. In the following subsections, we shall briefly discuss 
some of the earlier work on adaptive mesh generation strategies. 

4.2.1. Adaptive strategies by Finite Difference Methods : 

For finite difference applications, several authors ( Rai and 
Anderson, 1982; Brackbill, 1982; Dwyer and Sanders, 1983a; 
Gnoffo, 1982; Nakahashi and Deiwert, 1985; Saltzman and 
Brackbill, 1982; and Thompson et al., 1985b) have employed the 
structured grid generation techniques based on curvilinear 
co-ordinate lines, in their adaptive schemes. In these methods, 
the grid point locations are determined with the goal of 
minimizing some measure of error such as truncation error. The 
grid points are moved along the curvilinear co-ordinate lines till 
the error measure or the gradient is equi- distributed. 

Some researchers (Anderson, 1983a, b; Rai and Anderson, 1981; 
Greenberg, 1985; Eiseman, 1983,1985) have proposed dynamically 
adaptive systems in which the grid points move in response to the 
development of the physical solution on the existing grid. Several 
problems in heat transfer, flame propagation and fluid mechanics 
have been considered in these works. The grid point distribution 
over the field readjusts dynamically to concentrate nodes in 
regions of larger solution variation as they develop, without 
reliance on prior knowledge of the location of such regions. The 
grid point distribution is also expected to retain a sufficient 
degree of smoothness and not lead to an excessively skewed grid, 
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else the truncation error will increase (Must in, 19/8; Thompson, 
1982a, 1984, 1985d) . 

For the time dependent problems, some authors (c,hiu ct al . 
1983; Dwyer et al., 1982; Dwyer, 1983b; Gnoflo, 1982) have adopted 
the strategy in which the the grid changes at each time step, in 
this approach the grid and the numerical solution have boon solved 
separately. The time derivatives are transformed so that no 
interpolation onto the new grid is necessary. Gome methods, 
however, change the grid at selected time steps and here 
interpolation must be used to transfer the values, from the old 
grid to the new grid. Both the approaches have been discussed in 
detail by Thompson et al. (1982b). The variable used for grid 
adaption may differ from problem to problem; for instance, the 
grid could be adapted according to the velocity gradient in a flow 
problem and the temperature gradient in a flame propagation study. 
Equi- distribution concepts based on variational methods have been 
widely used in the recent past for grid adaption and those have 
been summarized in many survey papers and books (Hawkon, 1987; 
Thompson et al . , 1982b, 1984 ,1985a-c; Sengupta et al. 1088). 

4.2.2. Adaptive strategies by the Finite Element Methods : 

For finite element solution techniques, well developed 
mathematical theories are available for estimating the spatial 
discretization error in elliptic/ parabolic problems. For time- 
dependent problems, space-time finite element schemes are also 
available (Oden and Demkowicz, 1986 b, c; Babuska and Rheinboldt, 
1982a; Davis and Flaherty, 1982; Gelinas et al, 1981 ). In the 

space- time approach, the mesh generation procedure depends upon 



the minimization of the L 2 -error of the computed solution at a 
particular time. However, this procedure is computationally 
expensive as it involves the solution of non-linear equations. 

Miller and Miller (1981 a, 1981b) have developed the moving 
finite element method in which the grid point locations are solved 
as additional variables by invoking the Galerkin formulation. In 
this approach, the residual is required to be orthogonal to all 
the basis functions of both the solution and the grid. Also, the 
relative motion as well as the minimum separation distance between 
grid points are controlled up to some extent. 

Dupont (1982) has studied the finite element methods for 
which the underlying function spaces change with time. The error 
estimates and their convergence properties have been well- studied 
for the smooth solutions of parabolic problems in one space 
dimension and have been applied to fluid mechanics problems. 
Further, the types of continuous and discontinuous mesh 
modifications have been examined and certain symmetric a priori 
error estimates have been provided. 

In the context of FE solutions, error estimates with certain 
desired properties have been developed for many classes of 
problems and these have been employed in the design of effective 
grid adaption procedures ( Babuska and Rheinboldt, 1978a, b; 1982a, 
Babuska and Gui, 1986; Bank and Weiser, 1985; Gago et al., 1983). 
A posteriori error estimates based on the residuals have been 
proposed by Babuska and Rheinboldt (1978 a, b) for linear 
elliptic problems and have been widely used by many authors in 
finite element analysis. But for non-linear problems whose exact 
solution is not known, it can not always be guaranteed that the 
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error is bounded by the residual in some appropriate chosen norm, 
although such a feature is normally observed from the point of 

view of computations. 

Several authors (Babuska and Dorr, 1981b; Bn busk a et al., 
1981c, 1984; Demkowicz et al . , 1989a; Gui and Babuska, 1986a,b,c; 
Oden et al., 1989; Szabo, 1986; Zienkiewicz et al . , 1989) have 

worked on h and h-p refinement techniques for the solution of 
linear elliptic boundary value problems. These refinement 
techniques are commonly used in finite element modelling, for 
enhancing the convergence properties of the solution and reducing 
the error in the solution. Besides, these adaptive procedures 
allow not only the final global energy norm error to be well- 
estimated , but in addition, give a nearly optimal mesh. The 
desired accuracy can always be obtained with suitable h-p 
refinements and the rate of convergence of the adaptive h-p 
version analysis has also been analysed. It has. been suggested 
that for high accuracy, the use of p version or h-p version 
adaptive analysis with a simple error estimator is most 
appropriate (Gui and Babuska, 1986 a, b, c; Zienkiewicz et al., 
1989) . Gui and Babuska (1986 a, b, c) have shown that in h-p 
refinement, it is possible to achieve an exponential rate of 
convergence to the exact solution for a large class of problems, 
even those including corner singularities. The p version usually 
shows a higher rate of convergence than the h version and its 
implementation in different adaptive situations are well 
documented. However, the practical implementation of an optimal 
h-p strategy is difficult and has not yet been achieved - 

Bieterman and Babuska (1982a, b) have studied a variation of 
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the classical method of lines (MOL) in which the space— grid is 
updated in an adaptive procedure. In their adaptive FEMOL 
algorithm, grids change discontinuosly in time as nodes are 
added or deleted. The primary goal of the adaptive FEMOL is to 
keep the space discretization error less than a user specified 
tolerance at all times. The a posteriori error estimator used here 
is similar to that proposed by Babuska and Rheinboldt (1978a) , for 
linear PDE ' s . The most important part of the grid construction is 
a pattern recognition procedure which has been used to determine 
the "shape" of the grid. The process is made efficient, by 
controlling the error and monitoring the ODE time step-size 
information. 

4.3. Problem Description : 

In this chapter, the error estimators developed by Bieterman 
and Babuska (1982 a, b) have been used to simulate one dimensional 
front propagation problems. Both h and h-p methods have been 
employed and the equi-distribution of a posteriori error has been 
achieved through the statistical procedure proposed by ( Carey and 
Oden, 1984a; Oden and Demkowicz, 1986a) . Detailed comparisons are 
presented regarding the grid adaption capabilities of the h and 
h-p version strategies in terms of error behaviour, computational 
time/space requirements etc. The results from these two strategies 
are also compared with those of an inexpensive grid generation 
procedure based on the grid speed approach. 

Two specific application problems have been considered here 
which display the feature of a sharp moving front in a stationary 
medium. The first is a linear test problem with a known source 
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equations (4.1a-c) is given as below. 

u(x,t) = 1/2 + 1/2 tanh [2/3 (x - Ct) ] (4.2) 
This solution represents a wave front with approximate width 
(2/3 x ) that moves in the positive x - direction with a speed C, as 
shown in Fig. 4.1. Apart from a thin moving front, it also 
includes decay and sharp transition regions in time and thus, has 
a form similar to those arising in flame or shock front 
propagation problems. Bieterman and Babuska ( 1982a, b) have 
implemented a variation of the classical method of lines (MOL) for 
the solution of the above type of problems in which the space grid 
is updated adaptively. 


4.3.2. Description of the Flame propagation problem : 

One- dimensional pre-mixed flame propagation in a stationary 
gaseous mixture is considered. The details of the formulation have 
already been explained in the Chapter 3 and here the governing 
equations are recalled, for the sake of convenience. The transient 
heat and species diffusion equations for the flame propagation are 
given in a normalized form as : 


d T 

a t 

8 P 

a t 


a 2 T 
a x 2 


Ip 

a x 2 


= f(p,T) 

= - f(p,T) 


t > 0, X € (0,1) 


(4.3) 


with the same boundary and initial conditions as mentioned in 


equations (3.1b-c) of Chapter 3. The above equations can be 
expressed in a more general form as follows: 



Fig. 4.1 TEST PROBLEM SOLUTION PROFILE AT DIFFERENT 
TIMES 
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(u ! } t ” ( a ! (x) ( u ,) x ) = (x, t) eQx (0,t f ], 

1 = 1 , NPDE 

(4.4) 

For simplicity, we shall denote the Test problem (i.e. equation 
4.2 ) as problem (T) and the flame propagation problem in the 

absence of flow (i.e. equation 4.4 ) as problem (F) in this 

chapter. 


4.3.3. Adaptive grid notation : 

Here, we briefly describe the notations employed to solve the 
evolution problem on a given time interval, say [t , t ], by the 
adaptive h and h-p refinement algorithms. For the sake of 
convenience, we consider the sub interval [t , t ] of the 
global time interval [t , t ] and present the notation used for 
all the necessary variables over this time step. Here, the 
subscript variable i in the quantity t takes values 0,1,2, ...,f 
where t is the starting time and t f is the final time of the 
evolution problem. 


Approximation Spaces and the Solution ; 

First, we recall the finite dimensional solution space 
S r,k (A) which contains the approximation to u (x,t) of the 

h i 

equation ( 4 . 4 ) for each 1 (1 £ 1 npde). It is evident from the 

earlier discussion that the linear degree approximation of finite 
element solution space, i.e., 


,k 


(A); r = 1 , and the vector k = ( 1 , 1 , ... 1 ) 


(4.5a) 


is to be considered in the h refinement adaptive grid algorithm. 
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For simplicity, we shall denote this finite dimensional solution 

space i.e. , 

S^’ k (A) as y (A) . (4.5b) 

In the time interval [t , t ] , we write 

/ ,m (A) (4.5c) 

for the finite dimensional solution space at the m' !l adaption 
at t = t , which contains the approximation to the solution 

u } (x,t) for t e [t , t j+ ] and for each i { i i -r nphk) . 

Further, the dimension of y 1,m ( A) is A' 1,ra which represents the 

n 

number of basis functions and is the number of elements in 

e 

the computational domain A 1,m . The basis functions are denoted by 


i ,m 

0 (X), ( j = i,..,.jv 1 ’ ra ) ( 4 . 5d) 

, , n 


and /’ m (A) 


M 1 > m 

n 



f i ,m 

i n \ 

l , ih 5-*""*-. I,(ii \ „ m j 



II 

U l,j(X^) 

t, . (K) 


V. 

(4.5c) 


where U^^ft] are the nodal values of the approximate solution 

i,m 

u 1(j (x,t). Usually, with each approximation space y 1,w ( A) we 
associate a real vector 


= [ */• 


m l.« 


i , m 


... h 

N l ’ m 


(4.5f) 


where each component (i < js A^’ m ) represents the mesh 

size for the element fi 1,m ^ ' 1 ■ m i.™ .1 - 


= ( x 


i , m 

j ’ 


x l ’ m ) 

3*1 ' 


of the 


V « ) * * J 

computational grid A 1,m . Adaptive grid procedures typically 
employ a sequence of finite dimensional 


spaces 
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’ 1 (A), y l ’ 2 (A) , 


,y 1 ’ M (A) , 


with computational space meshes 


A 1 ’ 1 , A 1 ’ 2 , 


, A i,M . 


1 H . 

Here, P ’ is the space in which the final approximation is 

sought, with coarser spaces y*’ 1 ,/’ 2 In adaptive 

processes, however, there are no auxiliary spaces and the 
relationship between spaces is more involved. The nested 
approximation of spaces is not assumed because it poses severe 
restriction and is usually not met if there is node removal. 

In the following discussion, we may refer to the vector U i,m 

i , m 

by its component U (x,t) and the vector 


N l ,m 

i > m / i j ui i n 

U 1 [t]- { u ltJ [t]| 


(4 . 5g) 


Further, we denote by 


j = i 


uj’ m (x,t), 0J’ m (x,t) and uj ,ra (x,t) 
the computed solution, the initial condition and the starting 
guess solution available at the m th adaption at the time instanr 
t for each i (i ^ l £ npde) . 

We also use the following notation for the h-p refinement 
method. Let p be the n th component of the vector p defined as: 


P = 


V V 


p , }. 

y> m J 


( 4 . 5h) 


Here, P represents the degree of finite element approximation 
considered in the p version refinement of the n th element. In the 
case of h refinement components of the vector p (= u ) have the 
value as 1. In an analogous way, we define the sequence of 
finite dimensional solution spaces as 

y l,1 'P(A), y 1,2, P(A), ,y i,M 'P(A). 
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For the h~p refinement scheme the approximate solution:* f.or the 
initial condition, starting guess solution and the computed 
solution at t = t for m th adaption are denoted by 

6 i 1 ,m, P(x,t) , u|’ m ’^(x,t) and uj’ m ’P(x,t) 

The following notations about the errors will be used quite often. 
g 1,m (t) denotes the a posteriori error which is obtained by 
summing up the local elemental error indicators 7? J (t) for the 
n th element fi 1,m of the computational mesh A , ’"'(t) at t - t, for 

n * 


each i 

(i - 

l 

s npde) . Further, 

E 1 

n 

■ m (t) 

denote: 

5 the error 

indicator 

for 

the 

n t h element n 1 ’ m , 

n 

0 

£ n £ 

N 

1 , u\ ^ 

r 

and n 1 ’ the 

mean and 

cr 1>m 

the 

standard deviation 

of 

the s 

et 

O f 

local error 

indicators E 1 ' 

n 

m (t) 

; (i i n s N i,m ) 
e 

of 

all PI 

)E ' 

s at 

the instant 

t = t . 

The 

above symbols defined 

for the 

h 

and 

h-p version 


methods have been summarized in Table 4.1. 

Interpolation between Approxim ation Spaces : 

In the description of adaptive grid procedure, we shall need 
the interpolation operator 

q i „ m + l 
' 1 , m 

which is a linear transformation 

from y i,m to y 1,m+1 such that 5 ,,m+1 fu 1 ’ m (x, t) 

\ 

is " close " to the approximation U I,m+1 ( x ,t). In the nested 
or equal case y 1,m £ y l - m+1 , it is natural to take as 

1 * m 

the identity operator ? (although this is not always the best 

choice). The interpolation operator y ! > m+1 j. s use ful for the 

1 , m 

interpolation of solution variable from the m th adaption to the 

(m+i) adaption at each time t = 1 1 . In a similar fashion, we 


( . ) 
1 ,m ' ' 




9 


Table 4.1. Schematic representation of variables used in the 
h and h-p version algorithm 
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also denote ?| +1 as a linear interpolation operator f rom the time 
t = t to t which will be useful for time march inq in the h-p 
version. The order of interpolation may vary according to the 
circumstance. The global smoothness which depends upon the 
interpolation operator + 1 or 9 \* J is not considered in 

the h refinement scheme. For example, in the finite element 
context, if polynomial interpolation is applied independently 
for each function uj’ m (x,t), the smoothness is measured in 
terms of the smoothness of the individual functions which are 
defined on each element of the mesh. 


4.3.4. Finite Element Formulation : 

The details of the finite element formulation have been 


explained in the Chapter 3. Here, we state the necessary equations 
in the FE formulation which are useful in this chapter. The time 
Galerkin approximation of the solution of the system of equations 
(4.4) can be written as follows. 


a_uj(.,t) v a u (x,t) a v . ? rn . 

'at > V ' a i ax 'ax 1 ^ <f 1 (x,t,U),v i > 

for all U Vj e S^ ,k : t e (0,t f ] , and i « i , . . . , n m k . 

{ 4 . 6a) 

with suitably chosen initial conditions. The above system of 
equations (4.6a) can be written in a matrix form as 


W{ ^ | + [ K i] = { F x ( u t t 3)| for 1=1 ,NPDE 

(4.6b) 

with the appropriate initial values. Employing the Crank-Nicolson 
time discretization scheme the matrix system for (4.6b) can be 
written as given below for each t « r t t 1 

1 l—i ' i J ‘ 
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[ A l] { u ! [t] } = { g !< u1 [t]) - 

The details of the various 
have been discussed in the 


+ | H I ~ 1 (U i “ 1 [t]j.; (isis npde ) 

(4.6c) 

quantities in the equations (4.5a,b, 
Chapter 3 . 


c) 


4.4. Scope of the present study : 

4.4.1. Objectives : 

The basic aim is to develop a computational procedure which 

will meet the following goals : 

(Gl) . To obtain sufficiently accurate finite element solution 
to nonlinear partial differential equations representing 
flame propagation, where the accuracy is measured in some 
norm associated with the problem. 

(G2 ) . To obtain optimal mesh distribution at each instant of time, 
which accurately captures the shape and location of the 
moving front. 

(G3 ) . To develop h and h-p refinement strategies which meet 
certain tolerance criteria on the total a posteriori error 
for the whole solution domain. 

(G4) . To develop a grid speed based adaptive strategy which is 
computationally inexpensive. 

(G5 ) . To present a detailed study on the behaviour of error 
measures for the h and h-p refinement strategies and their 
application to Test problem and the Flame problem. 

(G6 ) . To compare the relative performance of the h and h-p methods 
with each other and also with that of the grid speed 
strategy for the two problems. 
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depend on controlling the spatial discretization error through the 
predicted solution for each time. Also, the optimal mesh 
distribution is computed based on the principle of 

equi-distribution of spatial discretization error with the help of 
elemental error indicators. The statistical procedure proposed by 
Carey and Oden (1984a), Carey (1988) and Oden and Demkowicz (1991) 
has been employed for achieving equidistribution of spatial 
discretization error. Since, the exact spatial discretization 
error is not known for the flame problem, it is replaced by the a 
posteriori error estimate proposed by Bieterman and Babuska 
(1982a, b) which involves the finite element residue over each 
subinterval of the mesh. But, for the Test problem in which the 
exact solution is available, the performance of the a posteriori 
error has been evaluated by the various error measures described 
in Chapter 3 . 

4.4.2. Mathematical formulation of the objectives : 

In this subsection, the objectives are expressed in the 
mathematical set up. Let us recall the a-priori, error estimate 
which is given by 

||( e (x,t) III = | { a ( x ) | x (e(x,t) ) 2 dx j- v t e (0,t f ) (4.7) 

The essential aim of the present work can be broadly classified 
into two parts. 

(a). Given a problem (T) with exact solution u(x,t) and a norm 
III III > design an efficient adaptive algorithm for constructing 
a finite element mesh 'A' such that 

(I! e ( x , t ) HI s TOL , V t € (0,t f ) (4.8) 

where U(x,t) is the finite element solution on the mesh A 
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and TOL is the given error tolerance. 

The a priori error j| e(x,t) J|| is replaced by a posteriori 
error fs'(t) since the present adaptive algorithm is based on 
&(t) . 

(b). Given a problem (F) and a norm |j| J|| , design an efficient 
adaptive algorithm for constructing a finite element mesh A 
satisfying equation (4.8) for each 1 ( 1 < 1 < npde ). 

The a priori error ||| e(x,t) ||| is substituted by a posteriori 
error & (t) since not only the present adaptive algorithm is 
based on S’(t) but also the exact solution of the system 
(4.4) is not known. 

The above aims (a) and (b) of the present work can be written 
in a simplified way which is as follows : 

The main aim is to control the a posteriori error 

g(t) s TOL , V t e (0,t ) (4.9) 

and the perform mesh modification in order to 
achieve at different time instants . 

Description for each time : 

The above aim has the following detailed description for each 

time interval; [t Q ,t i ] ; C t i ,t 2 ]; C t 2 ' t 3 ]; [t f _i ,t f ] of the 

total time interval [t ,t ]. 

o f 

The adaptive algorithm seeks to construct a mesh 

A i,M with mesh size vector ’ M 
and corresponding discrete FE solution 

u|’ M,p (x,t); ( i < 1 £ npde) 


such that 
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s' (t) =s TOL , for each t e (t ( ,t ( + i ] 
usually by constructing a sequence of meshes 


A 1 ’ 1 , A 1 ’ 2 , A 1 ' 3 , 


of mesh size vector 

1 1 T^i.2 1 3 ^ 

with corresponding solutions 


(4.10) 


(4.lla) 


(4.11b) 


l, 1 ' P (X,t), U^'PfX.t), U i ' P (X, t) ; 


U i 1 ,m ' p (x,t) e y 1 ,m,p (A) 


The mesh A 1,m+1 j. s constructed from 


U 1 ,m ’P(x,t) e y ,1,m,p (A) 


(4.11c) 


(4 . lid) 


by equidistribution of the contributions of different elemental 


error indicators 


E 1,m (t) ; (l s n s N i,m ) 

n cs 


to the a posteriori error 


r ,m (t) , t € [t . t 1 

v L i l+i - 


Here, the solution U 1 ,m 'P(x,t) is obtained by solving the semi 
discrete problem 


A r u;-p ( x,t) 


- v > + < a ,fe 1 


g_u 1 ’ m ’ p ( x , t ) a v 1 ’ m 
dx ’ dx 1 ^ 


= < f t ( x , t , U 1 ' m ’ P ( x , t ) ) , V | ' m " P ( x , t ) > 
V o; — p (x,t), v;—P(x,t) C y'--P(t); (x,t) < n X ( 


(1 a i s HPDE ) # and n = i .... ,h 

(4. lie) 

with appropriate boundary and initial conditions. Hero, all the 
components of the vector p may take different integral values 
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depending on the refinement scheme. Clearly, our basic goal has 
two ingredients : Firstly, we want the adaptive algorithm to be 
reliable in the sense that the error control (4.9) is guaranteed 
by the construction. Secondly, we want the algorithm to be 
efficient in the sense that, the constructed mesh A is nowhere 
over refined. In order to achieve this, a priori estimate has been 
used in a certain sense directly for the problem (T). 

4.5. h-version adaptive strategy : 

In this section, the h - version grid refinement algorithm is 
presented. Various steps involved in the finite element solution 
and mesh modification are described in detail. 

4.5.1. Description of the algorithm : 

The algorithm for the time dependent evolution of the solution 
and the mesh distribution is described here. At the beginning of 
each time step, the initial data involves the specification of (i) 
initial condition (i.e., solution at previous time instant ), (ii) 
initial mesh, and (iii) starting guess solution for the non-linear 
source term. After each adaption within a time step, the mesh 
distribution and the starting guess solution change, while the 
initial condition requires interpolation to the new mesh. Inner 
iterations are performed on each adapted mesh to ensure 
convergence of the non-linear finite element solver. On the basis 
of such finite element solutions a posteriori error is estimated 
with the help of local elemental error indicators and decisions 
for mesh modification are taken. If the total a posteriori error 
over the whole solution domain is less than a prescribed 
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tolerance, the grid adaption procedure if. halted and the solution 
at this stage is accepted as the computed solution for the updated 
time level. The sequence of steps performed for each adaption 
cycle are enumerated below. These have also been shown in the form 
of a flow chart in Fig. 4.2. 

Step 1 : Calculate the finite element solution at t = t. : 

Calculate the FE solution at t - t for m ,h adaption on 

the mesh A l,m with given initial condition and the 
starting guess solution t = t (initially with i i and mi) 
Step 2 : Calculate the a posteriori error estimate : 

Calculate the a posteriori error estimate f l,m at t. = t 

i 

for m th adaption with the computed solution U l ' ra (x,t). 
Step 3 : Check for the application of adaptive grid : 

1 , m 

if 8 (t) £ TOL then stop the adaptions, qo to st ep 5. 

1 , m 

if 8 (t) > TOL then go to step (4) for adaptive grid. 

Step 4: (a). Implementation of adaptive h-ref inement. .scheme 

(i). For for any fi*' m e A l ’ m , if hinioi. < R*‘ m (t) < n t>a 
then RETAIN the element as it. is. 

(ii) . For any pair of elements, (Q 1,ro ,£2 1,m ) if 

n n + i 

E n ’ (t) A E ^j(t) * mxmtol, then remove the common 
node between these elements. 

(iii) .For any Q 1,m e. A 1,m , if E 1,m (t) > n 1,m 

n n 

then REFINE the element according to statistical 

procedure. 

(b). Adaptive grid cycle : 

Generate initial data on the new mesh A 1 ’ ra> 1 
according to STEP 4a, increase the adaption counter 



Start 


, i 




Set m = 1 




1 



Fig. 4.2. Flow chart for h-version algorithm 
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(m) by one and go to step 1. 

Step 5 : Time marching : 

if t < t then set t - t + At; construct and 

f 1 + 1 i 

generate the initial data on the new mesh A 1 * 1 ’ 1 at 
t = t , reset counter (m) to 1 and go to step 1. 
if t = t then stop the algorithm. 

4 . 5 . 2 . Important features of the algorithm : 

Finite element solution : 

The finite element solution procedure requires initial data 
such as initial condition and starting guess solution at each 
adaption of a t ime step. In the adaption c y e .1 o , the initial 
condition fi 1,m (x,t) is obtained for the m' h adaption by 
interpolating the computed solution at the previous time level. 
For evaluating the non linear source term on the right hand side 
of the equations (4.5a) at t = t ( an starting guess solution is 
needed. We choose the initial condition 0 i ' ,, "(x, t) e y l m (A) as 
the starting guess solution 

| = u[’ m (x,t) e y l ' m (A) at t ® tj 
for the initial mesh adaption. For the subsequent adaptions (m >i) 
the freshly computed solution is used as the starting guess 
solution. 

With the necessary initial data, the finite element solution 
for the matrix system of equations ( 4.5c) is obtained at each 
t = t for each adaption according to the solution procedure 
described in section (3.3.3) of Chapter 3. Iterative schemes are 
required for solving the non-linear FE equations and execution of 
the adaptive mesh refinement, so that fully converged solutions 
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are not obtained in the FE solution procedure during intermediate 
stages of grid refinement. 


A posteriori error calculation : 

The a posteriori error estimate at t = t for the m th adaption 
of the system of equations (4.4) is given by 


a N 


1 , m 


1/2 


i , m 

s (t) = 


Z 1 , m 

E n <*> 


n = 1 


r N 


i , m 


NPDE 


Y _ ] Y * |71 ^’- (t)|2 


1 , n 


n =1 ^ 1 = 1 


1/2 


(4.12a) 

where, E*’ m (t) is the local error indicator including all partial 
differential equations on the n th element (= Q 1,m ) and V ,m (t) 
is the local error indicator for i th partial differential equation 
(4.4) alone. The local error indicator is given by the expression 


<’“(t) r 

1 »n 


if a 1 = 0 


i , m 
X 

r n+ 1 


12 a 


l 


I r i,m (x, t) I* dx if a" ± 0 

1 1 , n 1 1 


i . « 


(l < 1 < NPDE) ; ( 1 £ n ^ ^ i,m ) 


(4.12b) 


Here, a* = a { (x‘ ,m +x^”)/2 ), (i s i s npde) and the function 

r 1 ’ m (x, t) is the residual of i th partial differential equation 
i r n 

(4.4) on the n th element ( neglecting the discontinuities at the 
nodal points ). The .global residual of the i th partial 
differential equation (4.4) at t = t ( is given by 
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r|’ m (x, t) = 

a 1 (x) lu x (x,t) - f, (X|t,U (x,t)) 

( i ^ i -r n i> i) k ) 

(4.12c) 

_* > "i . t h 

where, U (x,t) is the computed solution of i PDE at t = t ( for 
m th adaption. The local error indicator is evaluated via eight 
point Gaussian quadrature in eq. (4.12b) for each element. 

Mesh refinement strategy : 

Decisions concerning the mesh modification at each iteration 
depend upon the space discretization error. If the a posteriori 
error evaluated for a particular adaption at a particular time 
step does not satisfy the tolerance criterion, then we identify 
the subregions of the space mesh where the error in the predicted 
solution is large and refine such zones. Further, after the 
refinement, the necessary initial data on the ref ined mesh are 
calculated with the help of interpolation schemes. The iterative 
computational procedure for obtaining the solution is repeated on 
each new mesh until a specified stopping criterion is met. The 
typical variations of the a posteriori error with adaptions and 
time step are illustrated schematically in Fig. 4.3. 

Statistical procedure and element divisions : 

The particular scheme used here assesses and refines the mesh 
according to a statistical procedure (Carey and Oden, 1984 a,b; 
Oden and Demkowicz, 1986c; Babuska and Rheinboldt, 1982a). One of 
the principle feature of this procedure is that a given element 


IF u i (x ' t) 



Time/adaption 
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with a large computed error, may be divided into several sub 
elements The number of subdivisions is proportional to the 
magnitude of the local elemental error. The refinement and 
derefinement procedures are carried out on the basis of the local 
error exceeding or falling below the mean value ot the local error 
indicators respectively. Elements with error less than the mean 
value are not refined. Besides, if the elemental error is a small 
fraction of the mean error for two adjacent elements, they are 
allowed to collapse into a single element by removing the common 
node. Such a statistical procedure involving simultaneous 
refinement and de-refinement at different parts of the solution 
domain is applied for each time step until the global a posteriori 
error meets the tolerance criterion. After a few iterations of 
this statistical distribution technique a well graded mesh results 
for each time step, while control on the total a posteriori error 
is provided by the proper selection of the tolerance limit TOI.. 

We shall now describe the details of the statistical procedure 
for the refinement of the mesh modification. 

(1) . Consider the local elemental error indicators E^ ,1B (t), 

(1 si s A'’* ’ m ) evaluated at the m th adaption of t = t $ , as 
the population to be ranked in the order of their magnitude. 
Taking the data set consisting of these error indicators as 
a random sample of given size drawn from a large population 
the mean and standard deviation are computed. 

(2) . The standard deviation cr i,ra is given by the relation 
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V 1/2 

^T( E*’ m (t) - a i,m ) 2 
j = i 


(A^ ,m - 1 ) 


(4.13) 


We classify the elements into various adaption sets based on 
the magnitude of the local error indicator E i,m (t) in 

n 

relation to the mean a i,m and the standard deviation cr l ’ m . In 
the following sets, the variable 'n' varies from i to ^ 1,m . 


9 = -( Q 1 ,m : mintol <E 1,m (t) < n" m 

I n n 


9 = { (Q i,m , Q 1 ’“) : E i,m (t) A E i,m (t) s mintol < 

2 [ v n ' n+l' n v ' n+1 v ' 


Q l ’ m : h 1 ’" 1 s 5 X 

n n j 


9 = ] £ 2 * 

4 n 


i , m 

* m 


E 1,m (t) ^ a. 1 ’ 1 " A h 1,m > 8 

n n 


- Q i,m : E i,m (t) — mintol A h i,m > 8 1 
n n J 


(4.14) 


where mintol is the prescribed tolerance for the removal of 
nodes in the domain A i,m and 8 is the prescribed minimum 
mesh size considered in the solution procedure. It is 
observed from the definition of these sets that all of them 
can never be empty. 

Define the "refinement intervals" on the distribution of 


elemental error indicators by : 


-r i, m i,m| x i> 

I 1 = n. ,n. +kcr , I 2 = n 

, I fc = £ a i,m + (t-i )k c 


+ k c r i,m ,a 1,m + 2 k <r i,m , 


a I,ra + (t-i )k <r l,m , t k a 


(4.15) 
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where k represents the fraction of standard deviation 
defining each refinement interval. 

The criteria for removing , retaining and refining the 
elements can be explained with the help of the sets defined 
above. 


(a) . REMOVE : 


A/* > 111 

For any n, 1 < n ^ ^ N 

the node which is common to 


if 


Q 


1 , in 


,n 


\ » m 
nil 


) 


both the elements 


c T then 
is removed. 


(b) . RETAIN : 

For any n, 1 =s 


£ N 


i , m 


if Q 


1 ,m 


G P V./ } 

1 


element fi 1,m is retained as such, 

n 

f /^i i > nil i * m \ / /-n i y m i ^ rn * ^ ts 

(Q , fi 9 ) nor ( Q , Q ) e 

v n-l ' n ' v n ' n+1 ' 2 

(c) REFINE : 

For any n, 1 < n :£ N I,m , if SI 1 '" e T 


then the 
provided neither 


then we calculate 


■i for which E 1,m e I. and the element is ref i nod / times 
n i 

by the statistical procedure. 


(d) MINIMUM GRID SPACE : 

For any n, 1 < n ^ N 1,n \ if Q 1,m e P then the element 

e nT? 

is not considered for refinement. 


( 6 ). Computational implementation : 

Now, we implement all the steps (a) , (b) , (c) and (d) as 

described above for all the elements in the domain A l,m . 
First, the unnecessary nodes on the mesh A 1 ’ ro are 
systematically removed applying step (a) . Some elements of 
mesh A 1 ,ra are refined according to step (c) and the degree 
of refinement is based on the statistical procedure as 
explained in the step ( 4 ) . in the computational experiments, 
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the degree of refinement (i.e., the number of divisions 
introduced ) is selected on the basis of the magnitude of 

t h 

*(n) for the n element which is defined as 

*(n) = (E*’ m - /i i>n )/ k o* 1 ,m (4.16) 

Besides, some elements on the mesh A 1 ’ m have been retained 
according to step (b) and step (d) . Note that the 
restrictions have been placed on the minimum spacing of grid 
points in order to avoid the collapse of the mesh. In the 
actual computations, the value of k has been taken as one 
and the value of the ^(n) has been restricted to small 
numbers (say 1 or 2 ) in the wave front region in order to 
avoid the dense refinement on the mesh A 1,m . 

In the statistical procedure, the mean values of the 
successive residuals or error distributions regress toward 
zero and simultaneously a natural gradation of the mesh is 
produced. The strategy gives equal weight-age to the error 
over small intervals with high gradients and to the error 
on a much larger element at low gradient in less sensitive 
regions. In practice, the scheme is not sensitive to the 
type of problem considered. 

"til. 

Construction of refined mesh ( A 1 ’ 1 ”* 1 ) data at (m+i) 
adaption of t = 

In the beginning of the transient problem ( t e ( t Q , t^) the 
initial condition on the initial mesh is known. This initial 

condition is chosen as the starting guess solution for solving 

the non-linear system of equations (4.5c) at t = t^ Let, the 

initial condition and the starting guess solution be denoted by 
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e, l,m (x,t) e 


y 1 • m (A) ; uj ’ (x, t) c //' 1 • m (A) at t t 


1 o 

(4.17a) 

For the finer meshes ( m > i ) , the computed finite 


element solution for the 


adaption U ] (x,t.) r 


(A) on 


I y m 


the computational mesh A 
newly refined mesh points of A 


t is .interpolated on the 


at t = 

l,m+ in the adaptive scheme to 


obtain the starting guess solution U 


1 , m + 1 


( x , t ) 


„ i , rn + 1 . . , 

J (A) . In 


other words. 


Uj’ (X,t) 


3 l ’ m+l U 1,m (x, t) 

1 , m 1 


(4.17b) 


where + 1 i s a linear interpolation operator from the m th 

i , m 

th 

adaption to the (m+1) adaption at the time instant t . In a 
similar fashion, the initial condition on the finer meshes (m > i) 
as can be obtained from 

0 i,m+1 (x,t) - + 1 0 ' ’ m ( x , t ) (4.17c) 

1 l » m 1 

respectively. In terms of adaption sets, on a particular element, 
the initial data can be retained as follows.. Let. an element 


Q 


1 , m+1 


( 1 < n ^ N 


l , m+1 


) of the "refined mesh A 


\ „ m + 1 


Then 


the following holds 


If 

V e Q 1,m+1 

n 

e A 1,m+1 , 

then 




II 

8* 

m 

Q l ; m e Q l 

n* 

,m for some 

n’ 

( t 5 n'S ;V 1,m ) . 

f* 

Also, 

<c e Q 1,m 

e T or 

T or T or 

T 

or T . 


n’ 

1 

2 3 

4 

s 


In reality, x e may or may not be a nodal point of A i,m . 

Thus, the initial condition can be written as 
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i ,m+ 1 

ft, (V/t) 


ft 


i , m 


(i,t) if x 6 n 1 ;" 1 € 3° or T or 9 or T 

n 1 2 3 5 


N 


i , m 


r i , m+ 1 

i , m 


( 1 s 1 < npde) 

( [t] (i) 


n * = 1 

if i 6 £ ^4 an< ^ ( 1 - 1 - npde) ( 4 . 17d) 

In a similar fashion, one can write the starting guess solution 
as follows : 


i , m+ 1 

u i (v,t) 


i ,m 


u 


(m,t) if me n i,m £ T or T or 9 or 9 

V ' ' n’ 12 3 5 


N 


I , m 


(l < 1 < npde) 


?; ,m+1 ( u;’ m ( £,t) ) =y u|’“ [t] & 

i ,m 1 / t 1 » n 1 y n 

n ’ = 1 


if x e Q i,m £ ? and (i s i < npde) 

n * 4 


(4. 17e) 


Here, •? 1,m+1 is the linear interpolation operator which is widely 

i , m 

used in the finite element interpolation theory. Thus the 
calculations in the eqns. ( 4.17d-e) are repeated for all the 
nodal values of £ e A i,m+1 and the initial condition and 
starting guess solutions for (m+l) th adaption are obtained. 


Time marching : 

In the h- refinement technique, the approximation space 
considered is a finite dimensional space of piece-wise linear 
degree functions. Hence, in time marching, the information 
available at the preceding time step is fully utilized for the 
generation of necessary initial data. Thus the finally computed 
solution is given by 
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U i,M (x f t) 6 / ,m (A) at t = t on A 1,m of [t ,t ) 

I ! I I 

(4 . 18a) 

which is taken as the approximation for the .initial condition 
0 1+1,1 (x,t) e y 1+1 ’ 1 (A) at t = t on A 1 * 1 ’ 1 of [t ,t. 


i+i 


[ t , t. ] 
L i i + i 1 


Further, we choose 


( 4 . 18b) 


A 1 ’ M s A 1+1 ’ 1 and /’ M (A) - y 1 * 1 ’ 1 (A). 

(4 . 18c) 

Together with the above set of equations (4.18a-c) which 
consists of necessary data for the computation o t. the given 
evolution problem at the next time interval, the h refinement 
algorithm is fully defined. 


As the present problem is nonlinear no analysis related to the 

choice of time step has been performed. However, very small time 

steps have been employed so that the accuracy of the solution is 

not affected by time discretization. In other words, we have 

assumed that the resulting transient ODE systems have been solved 
exactly. 


4.5.3. Computational procedure 

A computer code incorporating the h-version adaption algorithm 
has been written and validated for the test problem as well as 
flame propagation. A finite element formulation has also been 
loped to solve initial-boundary value problems for vector 
Y s of partial differential equations in one space dimension 
and time. The solution algorithm employs Crank Nicolson method for 
cretization with a piece-wise linear degree polynomial 
n for spatial interpolation. For the purpose of grid 
nd testing the solution convergence, we have assumed 
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that there is no time discretization error over the whole time 
span of the problem. 

At each instant of time, the finite element discretization 
leads to the solution of a nonlinear system of algebraic 
equations. In the present algorithm, a simple local iterative 
(Picard) method has been used to solve the matrix system of semi 
discrete equations (4.4). It is to be noted that, the fully 
converged solution on the grid is not at all computed prior to 
refinement. However, as soon as inner iterations satisfy a simple 
form of convergence on the given mesh, the grid refinement is 
performed. The convergence of the iterative method did not pose 
any problem, as we have considered very small time steps with the 
tolerance value for convergence set equal to 0.0000001. 

Following the standard finite element procedure, the element 
contributions for the matrix system of ordinary differential 
equations have been generated and used to calculate the 
corresponding global matrices. In the Crank-Nicolson time 
discretization, small time steps have been chosen in order to 
obtain accurate as well as good estimates of space discretization 
error. - In practice, however, one wishes to take as large a time 
step as possible to cut down the computational expense. But, 
larger time steps, can introduce some unwanted numerical 
oscillations into the solution, in addition to decreasing the 
solution accuracy. For the test problem, we have chosen uniform 
time steps of 0.000001 throughout the time evolution process and 
for the flame propagation problems two different time steps 
namely, t = 0. 00001 ( before ignition ) and t = 0.000005 (after 
ignition) have been used. The initial condition for each time 
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step, has been expressed by its pioco-w i so linear degree 
approximation. For time t = 0, the continuous initial data i$ 
converted into the approximate initial vector by rer.tr iction to 
the nodes Xj of the initial mesh. Besides, for each time instant, 
the solution on the preceding mesh is interpolated suitably to 
obtain an approximation to the initial condition on the newly 
introduced nodes of the refined mesh. 

In the test problem as well as the flame problem the following 
tolerance parameters have been considered through out the problem. 
For the de-refinement step in the adaptive procedure, local 
tolerance parameter ( REMOVE TOL ) has been chosen as 0.0000001 
for the removal of any node common to two adjacent elements, and 
the local tolerance parameter ( RETAIN TOL ) has been taken 
0.000001 for retaining the elements . Besides, the global a 

posteriori error tolerance ( GLOBAL TOL ) has been set equal to 

0.075, which is used for accepting the converged solution of a 
time step. In the refinement step of each adaption, a combination 
of one/ three/ seven new nodes have been introduced in various 
elements according to the statistical procedure for adaptive 

refinement during the ignition phase or start-up transient 
propagation phase. But, in the wave front / flame front region, 
the refinement is confined to one or three new nodes in any 

element. This restriction is considered in order to avoid over 
refinement. 

The error estimators and the nonlinear source term in the 
flame propagation problem have been computed by eight point 
numerical quadrature. The decisions concerning the mesh 
modification are made adaptively and automatically in the 
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programme at various stages of the operations. All computations 
have been carried out in double precision arithmetic in order to 
avoid the round off errors. 

4.6. hp-version adaptive strategy : 

4.6.1. Description of the algorithm : 

In the h-p method, the p -version and the h -version 
approaches are used alternatively. In a routine way the FE 
solution is computed at a given time and the a posteriori error is 
calculated. If the tolerance criterion on a posteriori error 
estimate is not satisfied at time t = for the m th adaption, 

then either p refinement or h-p refinement is performed. The 
algorithm identifies regions where the errors are locally large 
using statistical procedure. The problem is first solved with 
linear degree finite element approximation at t = t on the 
initial mesh. The de-refinement procedure is applied during the 
initial adaptions for various time intervals. The degree of finite 
element approximation is increased locally in the regions of 
large error and the computation is carried out till the a 
posteriori error tolerance criterion is met. But in case the error 
remains above the tolerance at t = t and at the same time the 
prescribed maximum degree of polynomial approximation has been 
already reached in some elements, then the mesh size h is reduced 
locally in the regions where the higher order FE approximation is 
still needed. It is done by introducing new grid points in these 
elements along with simultaneous replacement of higher order 
approximation by piece— wise linear degree approximation. 

After the construction of the refined mesh the problem is 
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solved with the linear degree approximation ot solution on the 
new mesh. The a posteriori error estimate at t t ( on the 
refined mesh is also calculated. If the a posteriori error 
criterion is not still satisfied, then the degree of finite 
element approximation is increased till the error is brought down 
the tolerance level. In this way p or h-p strategy is applied 
judiciously to obtain the adaptive mesh as well as the FE 
solution. 

As we advance in time from t^o t j + j the approximate solutions 
for the initial condition and the starting starting guess solution 
are both chosen from the finite dimensional solution space of 
linear degree approximation in each element over the initial mesh 
at t = t . If the converged finite element solution on the mesh 

I H 

A ’ at t = t (which satisfies the a posteriori error criterion), 

contains higher degree approximation in some regions, then these 

regions are locally refined and the corresponding FK solutions are 

converted into equivalent linear degree approximation. This 

procedure is carried out in a way, as it has been explained for 

the h-p refinement. The interpolated computed solution on the 

refined mesh or the actual computed solution on the A 1,M at t. " t ( 

is chosen as the initial condition at t = t the initial mesh 

1 + 1 

i + 1 i 

(A )* In the following, various steps involved in the FE 

solution and mesh modification are described in detail. The flow 
chart for the h-p algorithm is shown in Fig. 4.4. 


Step 1 : Calculate the finite element solution at t = t : 


Calculate the FE solution at t = t 

i 

the mesh A I,m with given initial 


i 

\ h » 

for m" adaption 
condition and 


on 

the 
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I Start I 

. T 

I Set m = 1 | 


1 



Fig . 4.4 


Flow chart for h-p 


version algorithm 
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starting guess solution (initially t i and m 
Step 2 : Calculate the a posteriori error estimate : 

Calculate the a posteriori error estimate t 1,m (t) with 
the computed solution U 1,m (x,t). 

Step 3 : Check for the application of adaptive grid : 

i , m 

if & (t) £ Tol then go to step 6 for advancing in time 

l , m 

if & (t) > Tol then go to step 4 lor the application 

of h-p refinement. 

Step 4 : Implementation of hp refinement, scheme : 

For each element ’ m is,, .< A l,m and do the following 

operations . 

(a) . For any Q 1,m e A i,m , if E l,m (t) < t ! '"‘ r hen 

n n 

go to step 5a for de-refinement procedure. 

(b) . For any Q i,ra e A 1,m , 

rt 

if E* ’ m ( t) > and Deg( S2 t,m ) a MAX DKG then 

11 n 

go to step 5b for the application ot p ret inement. 

(C). For any Q 1,m e A 1,m . 

n 

if E^ ,m (t) > n l ' m and Deg( u‘* m ) MAXlfKG then 
go to step 5c for conversion of p to h refinement 

Step 5 :(a). Application of de-refinement scheme. 

The two operations RETAIN and REMOVE of nodes is 
done according to step (4) of h-ref inement algorithm. 
After derefinement, go to step 6c. 

(b) . Application of p-version : 

For for any n‘' m (t) e A 1>m , if Deg< Q 1>m ) < MAXDEG 

ft 

then increase Deg (Q^’ m ) by l and go to step 6a 

(c) . p version to h version conversion : 

lil- If Deg (SI 1 ’*) a i and E !,w (t) < a 1, 

n n 


then 
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retain the element configuration and consider the 
piece-wise linear degree approximation in Q i,m 

n 

lAi-I * If De 9 (n i,m ) = 1 and E i,m (t) ^ n. 1 ’ 1 ” then 

n n 

apply h-refinement by the statistical procedure, 
(iii). If Deg (Q i,m ) > maxdeg and E i,m ft) ^ n i,m 

n n 

then refine the element into piece-wise linear 
elements by the p to h version conversion, 
according to step 6b. 

Step 6 : Construct the refined mesh A 1,m data : 

(a) , p version : 

Increase the adaption counter (m) by 1. 

Incorporate the additional degrees of freedom and 
construct the initial data on A i,m . Then go to 

step 1 . 

(b) . p version to h version : 

Increase the adaption counter (m) by 1. 
generate refined mesh A i,m+1 with piece wise 
linear FE approximation by p to h conversion, 
and go to step 1 for next adaption. 

(c) . De-refinement : 

Increase the adaption counter (m) by 1. 

construct the initial data on A ’ . Then go to 

step 1 . 

Step 7 : Time marching : 

(a) . If t ( < t f , 

For each Q i,m (t) e A l, \ if Deg( n‘ ,m ) > 1 

n 

then refine J2 l,m by p to h conversion and obtain 
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the piece-wise linear degree approximation in the 
refined elements of Q 1 ’ . Generate the initial. 

n 

data on the refined mesh A 1 " 1 ’ 1 and go to step l. 

(b) If t = t then stop the algorithm. 

4.6.2. Important features of the algorithm : 

Hierarchical forms of higher order elements : 

The hierarchical basis functions considered in the present 
study, simply represent additive refinements of higher order over 
the linear degree FE approximation. Let us consider the generation 
of these basis functions for a general element as illustrated in 
Fig. 4.5. 

The approximation over this element cannot be improved upon as 
the identification of the connecting nodes between elements 
guarantee C° continuity. However, we can achieve the hierarchical 
format over this expansion by preserving the required C° 
continuity of the linear approximation between elements. This 
condition can be satisfied by choosing basis functions which 
vanish at the end points of the interval. A particularly 
attractive class of polynomial functions which preserve C° 


continuity at the end points 

are 

related 

tothe 

integrals of 

Legendre polynomials p (C) over 

p 

the 

range - 

1 s C 

s l. Here we 

define the Legendre polynomial of 

degree p by 



V (C) = - — 

p 

d 

(C 2 -i) p 



P p ! 2 p 

dC P 



(4 . 19a) 

and integrating 2 P p (£) for 

P 

p = 

1,2, 3, 4, 

from 

-1 to G in 

turn gives 3 ,3 , 3 , and $ 

etc 

. , which 

. are 

used as basis 

functions of quadratic, cubic, 

fourth 

and 

fifth order 
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Fig 4.5. Hierarchical element shape functions of 
nearly orthogonal form 
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approximation respectively . The structure ol the elemental 
matrices arising from these hierarchical basis functions i s 
simple, due to the orthogonality property of; the Legendre 

functions over the interval -1 s < l. 

By means of such hierarchic nesting of the basis functions, one 
is able to utilize the element calculations of previous adaptions 
repeatedly up to some extent. The key feature to be observed here 
is that the addition of higher order terms does not alter the 
entries corresponding to the lower order basis functions. This 
implies that a large block sub matrix from the previous level 
remains unchanged. Using such hierarchical families of C° 
elements in one space dimension, the approximation for the 

unknown solution variable 11(C) on a typical element Q for any 

partial differential equation can be written as follows : 

V- (C) = 11^(0 + u?(C) + Z £.(C) 

J - :» 

c € Q = [-1,1] 

(4.19b) 

where $ and $ are the linear basis functions and .6 ; i > 2 

12 ) 

are correspond to higher degree Legendre polynomials and 
are values of the unknown variable corresponding to the nodes i 
and 2 respectively. The additional degrees of freedom 1 

(1 - J - p) enter as tangential derivatives at the mid side 
nodes, the nodal co-ordinate data on the end points of all 
elements is remain. 


Finite element solution : 

In the adaptive grid cycle, with the given initial condition 


0. 1 , tn | U 
1 


1 I in j, XI 


(A) at t 


t . 

i 


(x,t) € y 
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and the starting guess solution 

u|’ ra ' p (x,t) e y i ’ m >P(A) at t = t [ 

I Hi 

on the mesh A ’ , the resulting non-linear system (4.4 ) is solved 
by the iterative scheme as described in Chapter 3. Here, the 
additional unknown variables introduced in the p-version have 
also been invoked in the FE formulation and solved together with 
the unknown variables u)’ 1 " [t] ( i s j s jV 1 ’”). These additional 

unknown variables which are present in the approximation of the 
starting guess solution are initially chosen to be zero ( the 
vector p = u ) . But for subsequent p version adaptions, these 
paramrters need not be zero. 


The a posteriori error evaluation : 

In the h-p refinement method, p refinement and p-h refinement 
are employed during different adaptions of a particular time step. 
During p refinement, the higher degree approximations have been 
considered for the evaluation of local error indicators -n i,m (t) . 

1 , n 

Thus the a posteriori error estimate £ i,m (t) is given by the 
relation (4.12a) in which the error indicator r? l,m (t) for the n th 
element is given by 


l,n 


12 


X 


i , m 
- n+ 1 


r |;: <>^>1 


dx 


i ,m 


if a n = 0 
i 

if a n ± 0 
i 


(4.20a) 

where p is the degree of approximation used in the n th element, 

n 

a" is defined as earlier. Further, the residue in the 


and 



A posteriorierror 
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Time/adaption 


Fig. 4.6a. Schematic representation of the evolution 
a posteriori error with adaptions/time in 
h-p version algorithm 


of 
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based on the mean and the standard deviation of local error 
indicators. We define the following adaption nets, in all the 
following sets the variable 'n' varies from i t « A ‘ ’ m 


T 

1 


9 

2 


f (n i ’ m ,Q i ’ ra ) : E 1,m (t) A E l, "(t) 

' n n+l' n n+ 1 

i , m 

Deg (Q ) and Deg (£2 t ) - 1 

V r\ n+ 1 

/ Q i>m : E 1>m (t) < a 1>m or h i>m < 5 

{ n n n 


MJNTOI. < n 


i , m 


P 


{ i » m 

fi 1,m : E i,m (t) >a 1,m A h 1,m > 5, Dog (£2 ) t 

n n n n 


T 


<P 


(t) ia. i,m Ah 1 ’ m >5, 1< Deg (£2 ) < haxii k r, 


= / n I,m : E 1 ’ 1 

In n 

r i , m 

= | : E^ ,m (t) 2 :a i,m Ah‘ ,m > 5, Deg(£2 n ) max u k i; 


(4.21) 


(2). For all elements, determine in which refinement: interval the 
elemental error resides and refine the element', accordingly 
either by p version or h version. 

In the computational experiments, the de-refinement has 
been done only for elements on which the liner degree of FE 
approximation has been considered. It was discussed earlier 
that the piece-wise linear degree FE approximation has been 
used m the present algorithm on the initial mesh as well as 
on the mesh obtained after the conversion of p to h version. 
On such meshes with linear degree approximation, the 
derefinement procedure of the h version algorithm is applied 
first as explained below. 
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De-refinement : 

(i) . REMOVE : 

For any n, i s if (Q i>a ,n i,m ) e 9 then the 

« n n+1 ' 1 

node which is common to both the elements is removed. 

(ii) . RETAIN : 

For any n, 1 £ n * N y ’ m , if n i,m e P then the element 

e n 2 

Q‘ ,m is retained as such provided neither ( fi*’“ ,0*’ m ) 
nor ( Q i,m , n 1,m ) € 9 . 

p version refinement 

i , m 

For an element Q e 9 on the mesh A 1,a and then 

n 4 

i , m 

increase Deg(Q ) by one and recompute the FE solution 

n 

on the refined mesh A? ,in+1 . 

Although, the refined mesh A i,m+1 = A 1,ra , but the degree 
of FE approximation may be different between the 
elements of A 1,m and the element of A 1,m+:L . 
hp version refinement : 

Here, we describe the construction of refined mesh 
A i > m + 1 from A i,m for piece-wise linear degree FE 
approximation of solution. Since some elements of the 
mesh A i,m are either linear or higher degree, we 
systematically apply refinement, and 
conversion of p to h version as discussed below. 

(i) . For n i,m e P we find the integer i such 

that that E i,m (t) e I. by the statistical 

n 

procedure as explained in h version algorithm and 
subdivide the element in proportion with i 

(ii) . If fi l,m € 9 then consider linear degree FE 



approximation and neglect the h Lghor degree 
approximation in fV’ m 

( iii ) . if fi 1,m e P then n* grid points are introduced 
in Q i,m in such a way that the higher degree 

n 

polynomial in q‘ ’ m is approximated uniformly by 
piece-wise linear polynomials on the newly formed 
sub elements of Q i,m 

n 

Now, we implement all the steps as described above for 
all the elements in the domain A 1 ’" 1 . 
conversion of p version to h version : 

The procedure for obtaining linear degree FE 
approximation from the higher degree p version over an 
element Q I,m and the conversion from p version to h 

n 

version is described below. Since the construction of 
new mesh is based on the local gradient of the computed 
higher degree solution, a non-uniform refinement of each 
higher order element is obtained after conversion. Here, 
only one substantive variable i is used. The- substantive 
variable corresponds to the the temperature variable in 

_i,m,p 

the flame propagation problem. Let U (x,t) denote the 
computed solution of temperature. Then 

_i , m, p i , ra , p 

Uj (x,t) e if may be the approximate p-version 


FE solution 

in Q 1 ' 

n 

m of 

A 1 ’ m 

. Then on A 1 ’ 

m we have 







^ » Wl f J) i 

n 

^ 1 i » nj 


i t m 

h 


Uj (X,t) = 


[t] 


(X, + P; ; : 

0 (x) 

> J ’ 


J=i j * — i 


x e A 1,m and i a l s ww: 


( 4 . 22 ) 



127 


In the above equation (4.22), the unknown parameters 
<xj’j[t] ) enter as additional degrees of 

freedom due to hierarchical p version approximation. The 
hierarchical basis functions ~5 (x) . y = 1 ,.., 

in this equation may be equal to quadratic/cubic/ fourth/ 
fifth order polynomials as described in the section 
(4.6.2). Further, note that the \p\' m (x) (is j s N i,m ) are 
linear basis functions and U* ,r "[t] (i^ j s N i,m ) are 

1 , j n, 

the corresponding solution values at the nodal points. 
For an fi 1,m £ A 1,m , having higher degree approximation 

n 

solution, can be written in terms of a linear basis as 
well as higher order hierarchical basis functions. 
Denoting the solution in element by 


S(x) 


_1 ,m,p 

U (X,t ), 

1 jf n * 


X € 



n 



X 


i ,m 
n+1 


(4.23a) 


We set out to find 


r M 


new grid points 


a, , a. 


a 4 , . . .a i 


• a M-i* a M 


(x 1,m s a s x l, “); (is J £ *) 

n j n+1 


and corresponding 
polynomials 
these subintervals 


piece-vise linear degree 
M ) defined on 


H 1 (x) ( o s j 


i ,m 


[x^’ m , aj, [a x , a 2 ], ' Ca ^' X n+i •* 

such that the polynomial ^(x), is uniformly 
approximated by piece-vise polynomials ?^(x) 

(o == j s M ) on nl’ m . That is 

s c V x e [ J 
a = x 1 ’ m and a 


*?(X) - Wj (X) I ^ c V x e [a jf a j + 1 3 and ° ~ J < M 


, i , m 


M* l 


(4 23b) 
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In order to obtain the new grid points satisfying the 
equation (4.23b), we proceed as follows : 

(i). First, construct the straight line £(x) passing 
through (x‘ ,m , 9 (x^’ ra ) ) and (x‘^" f •' ' ( x ’ " ) ) given 

by 


£(x) = 9(x ’") + (x - x ,m ) 


S(x I>m ) 

V n+1 ' 

- ??( 

+ - 
3 

1 

I , m 
X 

n 


X € Q 


(4 . 23c) 


(ii). Determine the stationary points of h'(x) - £(x) for 

x e n 

n 

find x* e Q 1,m -such that 

n 

( 5(x*) - £(x‘) ) = 0 = Vf(x') ( say) 

The resulting equation has the form 




1 »m 3 ~t , , 

“i ,ySZ V) 1 *’ 


!,'(x’’" ) - K'(x' ’ ”) 

v nfl ' n ' 


J’ =1 


(4.23d) 


Note that required grid points are nothing but the zeros 

of the function W(x) « 0. Since !?(x) is a 

polynomial W(x) is also a polynomial, with real 

coefficients and by Rolle's theorem there exists at least 
one real root of W(x*) = 0 which lies between x i,m and 
x n ;“ respectively. if n‘ ’ m « we shall restrict 

ourselves to only n* roots of W (x* ) = o depending on 
the degree of finite element approximation used in Q 1,m 

rx 

as given in following table (4.2). 



129 


* 

n = 

i 

1 

2 

2 

Deg ( Q i,m ) = 

n 

2 

3 

4 

5 


Table 4.2a Selection of gr 1 points in an element of 
p to h conversi 

$ 

(iii).The choice of these n grid points is done as 

follows : Let y j , y 2 , ' Y n’-i ,Y n> (n ’ ~ 1) be the 

calculated roots of W(x*) which lie between 

x and x of £2 . 

n n + 1 n 

* * 

(a) . If n' = n : all of these roots are taken as n 

grid points in the element fi^’ m 

(b) . If n.' > n* : Compute | - ^fy^l f° r 

* 

j = n' and choose the first n 

according to the values of |^(Yj) - ^(Yj) I i n 

descending order. 

(c) . If n' < n* : All the calculations in the step 2 

of mesh refinement strategy are repeated over 
the newly formed subintervals 


yj, [y x , y 2 ],.--,[y n /_ 1 Y n /3» lY n ' x n+i ] 

Now collect all the new roots and select the 


required roots (n. ) according to step (b). 


(iv). If y. 


y 2 - 


V * y * (n* z: l) be the required 

' * n - i fJ n 


i , m 


and 


, m 
n + 1 


(n * ) grid points which lie between 
of Q 1,m construct the piece— wise linear 


interpolation polynomials 
defined on these subintervals 


(x) ( o s j £ n ) 


(V). 


[y 0 ty„v V>' [ V y „‘.i 

•t m i y nt 

where y Q = x n ’ and Y n * +1 = x n + 1 * 

These piece-vise linear polynomials ?yx) (o £ J 


n ) 
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represent polynomial S’(x), x € ' closely. These 

steps described above for p-h conversion, have 
been shown schematically in Fig. 4.6b. 

The important features of the p to h conversion procedure can 
be summarized as follows. In order to derive maximum benefit of 
the earlier p version calculations, a gradient based non-uniform h 
refinement of the highest order elements has been implemented. 
Since the gradient is the driving force in the present applied 
problem, the region of large error in the high gradient region of 
the physical problem is taken care of. Further, the conversion of 
p to h version offers gain in the distribution of grid points. At 
the same time, the degree of FE approximation has been converted 
to piece wise linear degree and then is a consequent increase in 
the dimension of associated solution space. 

Construction of initial condition on the refi ned mesh at t •• t : 

_i , m,]p - 

Let 0 <x,t) e be the approximate FE solution of 

p-version on A 1 ’" 1 . Then the solution U 1,m, ^(x,t) can be expressed 
as a linear combination of linear and hierarchical basis 
functions according to the equation (4.22) with the coefficients 
u| ,r "[t] (i=s j < and a 1,m [t] (i sys //} ,m ) which 

correspond to the solution values at the nodal points and the 
additional freedom. 

Initial condition ; p refinement 

It is observed noted that the meshes Q 1,ro and Q 1 ■ m *" 1 are 
equal but the degree of approximation may vary in certain elements 



Element (e) 



p- version 
solution 


Converted 

h-solution 


Grid points of 
old mesh 


# New grid points 



Fig. 4.6b APPROXIMATION OF COMPUTED FE SOLUTION 
DURING p TO h CONVERSION IN A HIGHER 
ORDER ELEMENT . 
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from m th to (m+i ) th adapt ion . The procedure for the h version 

algorithm is implemented for obtaining the initial condition for 
nodal variables of the approximate solution on the refined mesh 
fl i,m+1 . Further, the additional unknown variables ( 1 [ t ] ' s ) 

which have been introduced due to p refinement, are taken to be 
zero for the approximation of initial condition on the refined 
mesh A 1,m + 1 . This means that we have considered only the linear 

i , m + i , u 

degree approximation of FE solution space if for the 

i , m * i , u 

approximation of initial condition 6 ) (x,t). Thus 

1 , m , U i , m+ 1 , U 

if = if V m - 1 , H 


and 


i , m+1 i , m+1 1 , m + 1 

dim ( y i,m+1,u ) = N +N = 

n h a 


where ’ m + 1 = o. 

h 


h-p version : 

Since the p version elements are converted to h vers. ion linear 
elements the recently computed solution on A 1 ’" 1 is interpolated 
over the refined mesh A 1,m+1 to get an starting guess solution 
at t = t ^ for (m+l) th adaption. Further, the purpose of, 
de-refinement, the starting guess solution is determined in a 
similar way as explained in the h version algorithm . 

i.e., for any £ e Q i,m+1 c A i,m+1 , we have 

n 

£ = «• € Q i | m c Q 1,m f or some n . ( I < n << jV 1,m } , 
n e 
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Thus, the starting guess solution can be written as 


i , m + 1 

U (V/t) = 

1 , n 


_ i , m 

u , , (£, t) if £ e Q i,m c p or 9 

1 , n n ’ 1 2 


( 1 ^ 1 ^ NPDE) 

N i,m 

e 


y 1 ’ m + 1 ( u 

i , m ' 

I;;,(*,t) ) 

= 





n » 

= i 


if 

£ e n 1 ’ m £ 

n * 

9 

3 

and ( i ^ i 

s npde) 

N l ’ m 

n 

i , m 


N i , m 
h 


XX > 

[t] ' . (<c~) 

+ 


( £ ) 

j = i 


j 

* = 1 


if 

£ e Q 1 ' m £ 

n ’ 

T 

4 

o r P and ( 

5 v 

1 < 1 < npde) 


(4.24) 

Here, we have made use of conversion of the higher degree to 
linear degree approximation, during the construction of the 
starting guess solution. 


Initial condition : 

The approximation for the initial condition 

i , m+ 1 ,11 i , m + 1 , U i , m + 1 , U i,m+l 

ft (x , t ) e if with dim if ~ N & 

at t = t for the m th adaption can be obtained in a similar 
fashion as described in the h refinement algorithm. 


Time marching : 

If the finite element computed solution 

i M p . “ 

U ’ ’ ( x , t ) e y 1,M,p at t = t ( on the 

the fine mesh A i,M belongs to finite dimensional solution space of 
linear degree in each qJ’ M (i * i * ) then we adopt the 
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same procedure for the determination of initial condition on the 
initial mesh A 1 + 1 ’ 1 at t = t j as in h version. Thus wo have, 

1+1 , 1 , U _1,M,U Ml.l.U 

G ~ (X,t) = U (X,t) € if 

On the other hand, if 

_I,M,p l.M.p ),M,p 

U (x,t) e if where if is a 

finite dimensional solution space of higher degree then for the 
elements of higher degree ’ M (i s i s jV' ,H ) of A 1,M the p to 
h version conversion is adopted and the solution is computed at 
the refined grid points by the interpolation of the solution 

_ I i M , p 1 , M , p 

U (x, t) € if 

Thus we obtain 

i + i , i , p i + i , i , p 

G ( x , t ) e if 

the initial condition at t = t on the coarse mesh Q 1 * 1 ’ 1 of 

i 

linear degree approximation. Naturally, the dimension of the 
solution space y 1+1 » x »P changes at various points whenever the 
computed solution with p version is obtained at the previous time 
interval. The initial condition for the refined mesh is chosen as 
the starting guess solution for the same mesh which is useful for 
the solution of the non-linear problem. 

4.6.3. Computational procedure : 

In this section, we summarize the important computational 
aspects of the adaptive h-p version procedure for the test problem 
and the flame propagation problem. 

The concept of hierarchical finite element basis functions has 
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been used for obtaining mixed order interpolation in the h-p 
version. Linear basis functions are first used on the elements of 
the initial mesh for each time step and then the order of 
polynomial approximation is increased (up to a maximum degree of 
five) in the regions of the mesh where the error is large. The 
identification of these regions is done by the statistical 
procedure. This leads to the introduction of higher order nodes 
and additional degrees of freedom for such new modes in some 
elements of the coarse mesh, although the number of elements 
remains unaltered. Thus, when the new degrees of freedom are 
introduced, the total number of unknown variables in the finite 
element approximation is increased. 

The maximum order of polynomial approximation has been set 
equal to five in the present work. Very high order polynomials are 
susceptible to unrealistic numerical oscillations and it was 
observed during computations that the quality of the solution 
deteriorated beyond p = 5. At each stage of adaption, the error 
indicators have been re-computed to test the error criterion. 

Since the refinement is made hierarchically, some of the 
matrices arising in the discretized form of differential equations 
(4.4) on the coarse mesh are preserved. This is a very useful 
feature which results in saving a lot of computational effort, 
while calculating the finite element solution. Further, the 
hierarchical degrees of freedom appear as perturbations on the 
original solution rather than its substitutes. So, the 
resultant matrix equations of the associated equation (4.4.) has a 
more diagonally dominant form than that obtainable in a direct 
approximation involving an identical number of non-hierarchic 
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degrees of freedom. This feature is very much useful. in 
accelerating the convergence of the iterative method employed for 
the resulting matrix equations. 

During the p version adaption, the computed initial solution 
on the coarse mesh with linear degree approximation has been used 
as the initial solution on the adapted mesh (with higher degrees 
of freedom) . Further, when the p version is converted into h 
version in the h-p refinement procedure, the new starting guess 
solution has been constructed with piece-wise linear degree of 
approximation on each element of the refined mesh. This is 
achieved by converting the local higher order approximation in any 
element to a piece-wise linear approximation on the newly formed 
sub-elements. In fact, the newly introduced nodes have been chosen 
in such a way that the piece-wise linear degree approximation 
closely represents the higher order approximation available on the 
master element. Now, the solution on the previous, mesh is 
interpolated to obtain the approximate solution at the newly 
introduced nodes of the refined mesh. The initial condition at a 
given time instant has also been calculated by a similar procedure 
whenever the final computed solution at previous time is obtained 
by the h-p refinement scheme. 

In the p refinement scheme, the degree of the approximation in 
the elements with large elemental local error indicators is 
increased by one and the solution is recomputed at the same time. 
Thus, the elements in the refinement region become hierarchic of a 
certain degree, whereas elements outside the refinement region 
would be hierarchic but of lower degree. The process is continued 
till the maximum prescribed degree of approximation is attained. 
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At this stage, if further adaption is needed for the satisfaction 
of the error criterion by the a posteriori error , we first 
classify the elemental error indicators according to 
classification of adaption sets as described in the section 
(4.6.2). Thus, the respective regions are identified and the 
refinement is performed. The tolerance used for the h-p refinement 
are the same as that of used for h-ref inement. 

4.7. Grid speed adaptive strategy : 

4.7.1. Description of grid speed technique : 

Several authors (Miller and Miller, 1981a, b; Dwyer and 
Sanders, 1983a; Larrouturou, 1986; Dervieux et al. , 1989; 
Benkhaldoun et al., 1988b; Ramos, 1985), have worked on moving 
grid methods and algorithms for the determination of the nodal 
velocities in the adaptive grid procedures. The concept of moving 
all the mesh points in the same direction at each time step has 
been worked out by Larrouturou and co-authors (Larrouturou, 1986; 
Dervieux et al., 1989; Benkhaldoun et. al, 1988b ) for the front 
propagation problems. In their work, the grid velocity C(t) is 
chosen in such a way that the integral of the solution variable on 
the computational domain remains constant. Ramos (1985) have 
calculated the flame speed using the maximum temperature gradient 
in the adaptive finite element method. 

In the following section, the development of a different and 
simple grid speed strategy has been discussed and the behaviour of 
various error measures during the application of this strategy to 
the test problem have been presented. In contrast to the previous 
works, here the. differential equation for temperature has been 
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used to calculate the flame speed after invoking the assumption of 
uniform front velocity. 


4.7.2. Description of the algorithm : 

In the test problem, the solution profile at each time 
instant consists of a narrow wave front region of large gradients 
and a wide zone where the solution is constant. As the wave front 
moves, the previously adapted grid lay out in the front region 
also has to move at the speed of the front in order to provide 
optimal spacing for later times. In the results presented earlier 
for h or h-p version schemes, error measures were checked for each 
time step and grid-adaptions were performed whenever necessary. 
This results in a large computational effort for the grid adaption 
process. A computationally inexpensive way of achieving grid 
movement would be to calculate the front speed for each instant of 
time ( from the computed solution profile ) and displace the grid 
points at the speed of the front. This idea is explained as 
follows. 

Consider a mesh with a fixed number of grid points, which are 
moving at each time step, i.e, the nodal locations x are 
functions of time. On a moving grid, the time derivative of any 
variable at a grid point can be expressed as : 


^ U(x (t) , t) 


SU dx 3U 

3x dt at 


where the nodal velocity ~ appears explicitly. In many 
situations, the solution profile does not change with respect to 
time and the same profile translates from one location to another. 
For such a case, in order to maintain the same relative position of 
a grid point within the front, the required grid velocity can be 
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calculated from the equation, 


0 


d 

dt 


U(x(t),t) 


5U dx au 
ax dt at 


The above equation implies that there is no change in the solution 
with respect to time, for a moving grid point. Now for the test 
problem, substituting for — from eqs. (4.1a), with a(x) = 1, 
and b (x) = 0, the grid velocity C(t) is obtained as : 


C(t) 


dx 

dt 


(f (x ,t) + ) / |E 

a* 2 ax 


(4.25a) 


The resulting new location for node j at t = t is qiven as 

i + 1 

follows : 


X 


j + 


(t ) = 

1 v i + 1 1 


X (t ) 
j V i } 


+ C(t. ) (t. 


i + 1 


V 


(4.25b) 


The important steps involved in the grid speed strategy are 
explained below algorithmically. The corresponding flow chart is 
shown in the Fig. 4.7. 


Step 1 


Step 2 


Step 3 


Calculate the finite element solution at t = t^: 
Calculate the FE solution at t = t on the mesh A 1,m 
with given initial condition and the starting guess 
solution t = t . 

Calculate the Grid Velocity : 

Calculate time derivative of computed solution at t = t. 
and spatial derivative of computed solution at t = 

...... au.su 

Calculate uniform grid speed ( C(t) ) = - ' dx. 

Generate New mesh at t = t. with grid speed . 

Calculate new mesh location of a node ’j’ at t = t. 
using x j + i (t i + i ) = + c ( fc i ) + i " 1 1 > * 
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P Start | 

Given x 1 (t ( ) , i 

' i 

generate mesh A 


Define initial data on A 



Calculate FE solution 
at t - t on A 1 


1 s 

Y or. 

V 

t ~ 

= t 


f 



Advance in time 
Set i • i+i 


Calculate grid speed C(t) 


Calculate new location lor node i 
x‘ +1 (t) = x 1 (t) + C t.) At 
and generate new mor.h A 1 * 1 


Calculate the initial data on A 1 * 1 
by interpolation of data from A 1 


Go to i 


Fig. 4.7. flow chart for Grid speed strategy 
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Step 4 : Generate initial data on the new mesh A i + 1 at t = t 

i+1 

generate the initial data on the A I + 1 by interpolation 
of solution on A to A 1 1 , increase the time counter 
by one and go to step 1. 

4.7.3. Computational procedure 

The computational implementation of the grid speed strategy 
for the test problem as well as the flame problem is same. 
According to this strategy, all the nodes (except those on 
boundary) have been provided the same velocity C(t) at a given 
time t. This velocity C(t) is calculated as the average of all the 
velocities at individual grid points. Computational experience 
shows that the average grid velocity in the wave front zone is 
very close to the exact grid speed ( = front speed ) for the test 
problem. It is also observed that as the grid points move with 
same velocity, collapse and stretching of elements occur near the 
boundary which affects the solutions. In order to avoid this, 
special treatment of grid movement is necessary near the domain 
boundaries . 

In the computational experiments, the boundary nodes of the 
domain have been locked and the interior nodes are allowed to move 
according to uniform grid speed. After a few time iterations, the 
nodes which are near to the boundary may shoot outside the 
boundary. So these which cross the boundary are reintroduced at 
the beginning of the mesh. Although, it is not necessary to 
reintroduce these nodes in the present problem, for more general 
cases it may be necessary. 

It is important to note that the oscillatory nature of the 



142 


solution profile may affect the evaluation of the gradient. 
Further, the second derivative of the computed solution has to be 
evaluated carefully as it has significant contribution in the 
evaluation of the grid speed in wave front region. Tlius, we have 
evaluated the first and second derivatives of the gradient by 
finite difference expressions which are second order accurate. The 
consideration of second derivative term in the evaluation of grid 
speed introduces some errors because of the choice of linear 
degree approximation employed in the finite element solution 
procedure. However, for second order accurate procedures, the 
effect of the linear degree FE approximation is negligibly small. 

4.8. Results and discussion : 

4.8.1. Application of adaptive strategies to Test problem : 

In this section, results of solution profiles, error 
estimates, and grid - distributions are presented for the flame 
propagation as well as the test problems. The test problem has 
been considered in order to analyze the effectiveness of the error 
measures used for grid - adaptions in both the h version and hp 
version methods. In fact, since the exact solution is available 
for the test problem, the true gradient errors and relative errors 
of the computed solutions can also be estimated. In view of the 
sensitivity of error estimates to the smoothness of the solution, 
a local smoothening procedure using piece-wise cubic spline 
interpolation has been attempted in h version for improving the 
estimates of the true gradient error. Detailed comparisons of 
the relative performance of h and hp version methods have been 
presented. Finally, applying the grid - 


velocity concept to 



h-version grid adaption, the construction of computationally 
inexpensive meshes has been attempted and the behaviour of a 
posteriori error on such meshes has been studied. 

(I). Uniform mesh results : 

Bieterman and Babuska, (1982a-b) have shown the equivalence 
of the various error measures (see, section 3.4, Chapter 3) for 
linear problems. Here, for validating the computational results 
of the present analysis comparisons of the computed and the exact 
solutions have been provided graphically. The characterstics of 
different error estimators have also been plotted. The results 
have been presented in the tabular form (table 4.3a-d) and the 
comparisons with the results of Bieterman and Babuska (1982a) have 
also been shown at closely approximate time values ( table 
4.3c-d). Further, the notations of the various quantities used in 
the subsequent tables have been described in the table 4.2b. 

The test problem was solved on a sequence of uniform meshes 
with 20, 40, 80 and 160 elements. The reason for considering the 
uniform meshes is that in adaptive grid problems, additional error 
may result if the movement of grid points is not synchronized with 
that of the front. Therefore, in order to analyze the nature of 
the error estimates without moving-grid effects, uniform meshes 
have been considered prior to adaptive grid solutions. 

(a). Comparison of exact and computed solution profiles : 

The exact and the calculated solution profiles obtained using 
uniform meshes for the test problem are shown in Figs. 4.8a-d, for 
different time instants. For all the meshes considered, the 
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Rem (Rm) 

Ref (Rf) 

Stf 

Adp 

N (Eles) 
Nodes (Nos) 
Dof 
Deg 

Apos . err 
com. grd 
Grd. err 
Eff.ind 
Non-lin-ter 
a-st 

comp sol 
Exact sol 


Number of elements removed 
Number of elements refined 
Nuber of elements stufed 
Adaption number 

Number of elements considered. 

Number of nodes in the domain 
Total number of degrees of freedom. 
Maximum degree of approximation used 
Aposteriori error 
Computational gradient 
Gradient true error 
Effectivity index 
Non-linear iterations 
a -strategy 
Computed solution 
Exact solution 
h-version 
hp-version 


Table 4.2b. Description of the variable 
used in subsequent tables. 
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Fia. 4.8 TEST PROBLEM {Uniform Mesh ) : SOLUTION PROFILE 


So ) ut ion 

Line type 

Cornp so In 

Exact so In 

— 


Adp 

Time 

Eles 

Nods 

1 

0 .000060 

20 

21 

2 

0 .001250 

20 

21 

3 

O .005000 

20 

21 

4 

0 015000 

20 

21 

5 

0 .025000 

20 

21 

6 

0 .040300 

20 

21 

7 

0 .055300 

j 20 

21 

8 

0 .070300 

! 20 

21 

9 

0 .085300 

, 20 

1 21 

10 

0 .097800 

20 

21 


Acfc 

Time 

Eles 

Dot 

1 

0 .000050 

40 

41 

2 

0 .001250 

40 

41 

3 

0 .005000 

40 

41 

4 

0 .015000 

40 

41 

5 

0 .025000 

40 

41 

6 

0 .040300 

40 

41 

7 

0 .055300 

40 

41 

8 
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40 

41 

9 
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40 

41 

10 
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40 

1 

41 








EXACT AND COMP SOLUTION EXACT AND COMP SOLUTION 



0.9 -// / ■ 

- ! 


0.5 H- > 


1.0 1.2 

N=160 d” 


0.5 4 


0.3 r 


o.i H 


/ 




Solution 

Line type 

Comp so In 
Exact so In, 

— 


Adp 

Time 

Eles 

Nods 

1 

0 .000050 

80 

81 

2 

0 .001250 

80 

81 

3 

0 .005000 

80 

81 

4 

0 .015000 

80 

81 

5 

0 .025000 

80 

81 

a 

0 .040300 

80 

61 

7 

0 .055300 

80 

81 

8 

0 .070300 

80 

81 

9 

0 .085300 

80 

81 

10 

0 .097800 

80 

81 


Adp 

Time 

Eles 

Dof 

1 

0 .000050 

160 

161 

2 

0 .001250 

160 

161 

3 

0 .005000 

160 

161 

4 

0 .015000 

160 

161 

5 

0 .025000 

160 

161 

6 

0 .040300 

160 

161 

7 

0 .055300 160 

161 
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0 .070300 

160 

161 

9 

0 .085300 160 : 

161 

10 

0 .097800 160 

161 
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Fin. 4.8 TEST PROBLEM (Uniform Mesh ) : SOLUTION PROFILE 
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Time 

Elements 

N = 20 

N = 40 

N = 80 

N = 160 

0.000050 

1.3396 

0.6598 

0.3323 

0.1662 

0.005150 

1.8644 

0.9412 

0.4700 

0.2348 

0.010250 

1.8785 

0.9419 

0.4708 

0.2352 

0.023600 

1.9847 

0.9423 

0.4708 

0.2352 

0.032050 

1.8534 

0.9496 

0.4708 

0.2352 

0.046550 

1.8171 

0.9479 

0.4708 

0.2352 

0.056050 

1.8126 

0.9417 

0.4708 

0.2352 

0.060050 

1.8990 

0.9444 

0.4708 

0.2352 

0.070050 

1.8990 

0.9444 

0.4708 

0.2352 

0.088050 

1.9540 

0.9396 

0.4708 

0.2352 

0.099800 

1.3472 

0.6760 

0.3335 

0.1664 


Table 4.3a. Aposteriori error predictions 


T mo 

Elements 

X XXtlfcl 

N = 20 

N = 40 

N = 80 

EB9 

0.000050 

1.4257 

0.6639 

0.3290 

0.1661 

0.005150 

1.9995 

0.9421 

0.4654 

0.2346 

0.010250 

1.9901 

0.9382 

0.4662 

0.2350 

0.023600 

1.5307 

0.8595 

0.4662 

0.2350 

0.032050 

1.1912 

0.9216 

0.4661 

0.2350 

0.046550 

1.4355 

0.8691 

0.4662 

0.2350 

0.056050 

1.7108 

0.8620 

0.4662 

0.2350 

0.060050 

2.0083 

0.9462 

0.4662 

0.2350 

0.070050 

2.0083 

0.9462 

0.4662 

0.2350 

0.088050 

1.2401 

0.9108 

0.4661 

0.2350 

0.099800 

1.4392 

0.6479 

0.3277 

0.1661 


Table 4.3b. Gradient of true error predictions 
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(i). Present work : 


Time 

Elements "I 

N - 20 

N = 40 

N = 80 

N = 160 

0.005150 

0.5478 

0.2581 

0.1275 

0.0643 

0.010250 

0.5450 

0.2569 

0.1277 

0.0644 

0.023600 

0.4192 

0.2354 

0.1277 

0.0644 

0.032050 

0.3262 

0.2524 

0.1277 

0.0644 

0.046550 

0.3931 

0.2380 

0.1277 

0.0644 

0.056050 

0.4685 

0.2361 

0.1277 

0.0644 

0.060050 

0.5500 

0.2591 

0.1277 

0.0644 

0.070050 

0.5500 

0.2591 

0.1277 

0.0644 

0.088050 

0.3396 

0.2494 

0.1277 

0.0644 

0.099800 

0.5268 

0.2372 

0.1200 

0.0608 


(li). Results of Babuska and Rheinboldt 


Time 

Elements | 

N - 20 

N = 40 

N = 80 

N * 160 

0.03200 

0.05600 

0.06000 

0.08800 

0.4325 

0.5697 

0.6437 

0.3473 

0.2634 

0.2431 

0.2818 

0.2673 

0.1308 

0.1309 

0.1309 

0.1309 

0.0649 

0.0649 

0.0649 

0.0649 


Table. 4.3c. Relative error predictions 


(i). Present work : 


Time 

1 Elements 

N - 20 

N - 40 

N a 80 

N * 160 

0.005150 

0.010250 

0.023600 

0.032050 

0.046550 

0.056050 

0.060050 

0.070050 

0.088050 

0.099800 

0.9324 

0.9440 

1.2966 

1.5560 

1.2658 

1.0595 

0.9456 

0.9456 

1.5758 

0.9361 

0.9991 

1.0038 

1.0963 

1.0303 

1.0906 

1.0924 

0 .9981 
0.9981 
1.0316 
1.0434 

1 .0099 

1.0100 
1.0100 
1.0101 
1.0100 
1.0100 
1.0100 
1.0100 
1.0101 
1.0177 

1.0006 

1.0007 

1.0007 

1.0007 

1.0007 

1.0007 

1.0007 

1.0007 

1.0007 

1.0023 

(ii). Results of Babuska and Rhei 

Lnboldt 


Time 

Elements 


N ~ 20 

N - 40 

N «• 80 

N = 160 

0.03200 

0.05600 

0. 06000 
0.08800 

1.3200 

0.7760 

0.8060 

1.6850 

0 .9950 
1.0760 
0.9030 
0.9610 

0 .9890 

0 .9890 
0.9888 
0.9870 

0.9960 

0.9960 

0.9960 

0.9960 


Table. 4 . 3d. Ef fect ivity ratio predictions 
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computed solutions correctly predict a nearly linear profile in 
the moving front region. However, the slopes of the computed 
solution profiles are not accurate for the coarse meshes ( N = 20 
and 40 ) because the steep gradients are not resolved accurately. 
Results of this study show that the shape of the computed solution 
for N > 40 coincides with that of the exact solution. In fact, 
for N > 40, the grid size is less than the approximate width of 
the front, which is an important requirement for accurate 
resolution of the gradients. The graphs of the solution profiles 
as well as those of the error estimators, indicate that the 
accumulation o.f spatial discretization error affects the timewise 
evolution also, in coarse meshes. The time discretization error 
itself, however, seems to be negligible for the small time steps 
used in the present computations. 

(b). Behaviour of a posteriori error, true gradient error, 
relative error and effectivity index : 

The results for the a posteriori errors, true gradient 
errors, relative errors and effectivity index on uniform meshes at 
different times are presented in the Figs. 4.9. and 4.10. The 
observed trends of the errors can be explained considering three 
separate stages of the front propagation, namely : the start— up 
transient, uniform front motion and the extinction transient. 
During start- up and extinction which occur near the domain 
boundaries, the a posteriori error is low (Fig. 4.9a). This may 
be attributed to the smaller contribution of source term in the 
governing equation and the smaller value of the front speed. In 
the uniform propagation stage, the a posteriori error undergoes 
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systematic oscillations due to the a L tomato under-prodiction/ 
over-prediction of the solution around the mean value, with 
respect to time. The average magnitude of the a post rriori error, 
however, remains constant at this stage. As the number of 
elements is increased, this mean value of error decreases linearly 
proportional to the step size h (or inversely proportional to the 
number of elements). Also, the oscillations in error with 
respect to time are eliminated. 

As regards the true gradient error ( Fig. 4.9b ) the 
following observations can be made. During the uniform 
propagation of the front, although the average magnitude of true 
gradient error is constant, the amplitude of timewise oscillations 
is large. This implies that the variation in the computed 
gradients with respect to time is more severe than that of the a 
posteriori error. However, as the gradients. arc? better 
represented through grid refinements, the oscillations in true 
gradient error vanish. An interesting trend seen is that the 
error in the computed gradients decreases in the shapp transition 
zones, although the wall gradients at the boundaries are not 
accurately represented for coarse meshes. Thus the minimum value 
of true gradient error occurs slightly away from the wall for the 
coarse meshes. It is also observed that for the linear test 
problem considered, the a posteriori error and true gradient error 
are almost of similar magnitude, although the true gradient error 
is slightly smaller for coarse meshes. On fine meshes, the two 
error estimates become exactly equal to each other. 

Fig. (4.10a ) shows the variation of relative error with 
respect to time. The behaviour of relative error is similar to 



Fia 49. TEST PROBLEM : TINE EVOLUTION OF ERRORS (UNIFORM MESH) 
(APOSTERIORI ERROR AND TRUE GRADIENT ERROR ) 
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Fig. 4.10. TEST PROBLEM : TIME EVOLUTION OF ERRORS (UNIFORM MESH 
(RELATIVE ERROR AND EFFECTIVITY RATIO) 
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that of the true gradient error. It is observed that the relative 
error decreases rapidly, as the mesh becomes finer . In Fig. 
(4.10b), the effectivity index which is the ratio of the a 
posteriori and the true gradient errors, is plotted against time. 
For the coarse meshes, the effectivity index fluctuates rapidly 
due to fluctuations in the a posteriori and true gradient errors. 
Obviously , the oscillations of the two error measures are not in 
phase. As the number of elements is increased, the oscillations 
in the two error estimates are eliminated and the effectivity 
index approaches unity. 

From the Figs. (4.9a-b, 4.10a-b), the following conclusions 
can be drawn. As the number of elements is increased, the 
computed solution approaches the exact solution and the a 
posteriori error closely approximates the true gradient error. 
Both of these errors go to zero linearly proportional to step 
size. Thus, for complex problems in which true gradient error can 
not be estimated, the a posteriori error may be used as an 
effective measure for grid adaption purposes, provided some 
necessary conditions on the PDE system are satisfied. Regarding 
the choice of the minimum number of uniform elements required for 
the test problem, it can be said that the step size h should be 
less than the front width (2(3) -1 . 

(II). Non-uniform mesh results ( h and hp refinement methods ) : 
(a). Solution profiles : 

The exact and the computed solution profiles for the h and hp 
refinement methods calculated on non-uniform meshes have been 
plotted in Figs. 4.11a and 4.11b. Also, the histories of grid 



154 


adaption at various time instants have been shown in t lie tables by 
the side of the figures. At each time .instant, the computed 
solution corresponds to the final adaptive mesh on which the 
solution meets the error criterion . In both the methods, the 
adaption is carried out by starting from a grid of 10 uniform 
elements initially. It can be observed from the figures that the 
computed solution is inaccurate at initial times (near the x = o 
boundary) and the differences gradually disappear as the mesh 
adapts at various times along with the front movement. The 
initial deviations are obviously due to the inadequate spatial 
resolution by the coarse mesh. For later times, however, both the 
shape of the front profile and the speed of movement match exactly 
with the exact solution for the h and hp refinement techniques. 

(b). Behaviour of a posteriori error, true gradient error, 
relative error and effectivity index : 

The variation of the a posteriori error with time in depicted 
in Fig. 4.12a. It can be seen from the graph 'that, the a 
posteriori errors for both the h and hp refinement methods are 
very high initially, indicating inadequate spatial resolution by 
the coarse initial mesh. But by the repeated application of the 
adaptive procedure, zones of large solution error are refined and 
the a posteriori error reduces rapidly and monotonical ly below the 
prescribed tolerance. As the solution front moves in time, 
however, the refined mesh obtained initially is not best suited. 
The a posteriori error therefore starts increasing until fresh 
mesh adaption begins when the error slightly exceeds the 
tolerance. After this adaption, the mesh is adjusted to the new 
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Fia. 4.11 TEST PROBLEM (h & hp version ) : SOLUTION PROFILE 
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location of the front and therefore, a posteriori error is 
considerably less. This process of grid adjustment need not be 
performed at every time step, since a posteriori error .in many 
cases may be below the tolerance level. Only when the error 
exceeds the tolerance ( in other words, the front has almost moved 
out of the earlier refined mesh zone ) , grid adaption is required. 
Due to this lag between front movement and grid adaption, the a 
posteriori error undergoes fluctuations with respect to time. A 
remarkable feature in the error analysis is that the new initial 
mesh at any time has been constructed with the best use of 
information available on the final mesh at the previous time. 
This reduces the number of adaptions in the total time interval. 

For both h and hp refinement methods, the adaption criterion 
works well to keep the error below the tolerance lor most of the 
time. In the case of hp method, the mean error is slightly larger 
than that in the h- version; however, the number of elements 
employed in the hp version is less. Further, it is observed that 
increasing the degree of finite element approximation in some 
elements locally, reduces the a posteriori error very slowly. 
Therefore, the number of adaption cycles per time step is more in 
the hp version in comparison to the h version. Close to 
extinction, the a posteriori error decreases for both the h and hp 
refinement methods, for reasons cited while discussing the results 
of the uniform meshes ( Fig. 4.9a). 

In Fig. 4.12b, the time evolutions of the true gradient error 
in the h and hp refinement methods are shown for the test problem. 
It is seen that during uniform front propagation, the mean values 
of the error are almost constant and equal for the h and hp 
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methods. However, some minute differences do exist between the 
true gradient error behaviours for the two methods . For instance , 
the amplitudes of the oscillations are smaller but the frequencies 
are higher , for the hp method . This is due the fact that the 
error decreases slowly as the order of interpolation is increased 
at the locations of large error; in fact, for higher order 
interpolations there is sometimes even an increase in the local 
error due to the pollution effect. Thus, more number of adaptions 
are often required in the hp version, which is shown in the form 
of higher oscillation frequency of the true gradient error. In 
general, the true gradient error for the hp method is larger than 
that for the h version. 

In the light of the discussions presented above, it is 
important to note that for initial times the hp method does not 
meet the error criterion even after reaching the maximum degree of 
interpolation in some elements. Therefore, essentially h version 
adaption is carried out in the hp version until the a posteriori 
error becomes less than the tolerance. Another important feature 
which is apparent from the comparison of Figs. 4.12a and 4.12b is 
that -the a posteriori and true gradient errors are close to each 
other, but exactly not equal, for the non-uniform mesh. In fact, 
the tolerance criterion serves as an upper limit for the a 
posteriori error while in the case of true gradient error, the 
mean value is closer to the tolerance. 

For non-uniform meshes, the relative error behaviour is 
similar to that of true gradient error as seen from Fig. 4.12c. 
The effectivity index ( Fig. 4.12d ) which is the ratio between 
the a posteriori and true gradient errors, decreases steeply 





0.010.01 0.03 0.05 0.07 0.09 O.il '0.0,1 0.0 d 0.03 0.05 0.0/ 0.00 0.11 


H 

on 

03 




159 


towards zero for initial tiroes and then increases slowly to reach 
a mean value around 0.7. Initially, the a posteriori error 
which is based on the residue is very large for the coarse mesh , 
while the true gradient error is small due to the small values of 
computed and exact gradients. Therefore, the effectivity index is 

i 

large. However, during the start-up period, the gradient at the 
wall increases (thereby resulting in larger true gradient error) , 
while the residue is brought down rapidly by grid refinement. Due 
to this, the effectivity index reduces almost to zero. In the 
front propagation stage, both the a posteriori and true gradient 
errors are of similar magnitude. So, the effectivity index 
increases. It is apparent from the Figs. 4.12a and 4.12b that the 
a posteriori error is more sensitive to grid refinement than the 
true gradient error. Thus, as the front moves, adjustment of the 
grid from time to time brings down the a posteriori error 
considerably but not the true gradient error. So, the effectivity 
index does not equal unity for non-uniform meshes. 

(C). The maximum and mean errors : 

In Figs. 4.13a and 4.13b, the variations of the maximum and 
the mean elemental error indicators with respect to time have been 
presented for the h and hp version methods. In the initial 
period, the mean errors are very high for both the refinement 
strategies, due to the coarse starting- mesh. Following the 
application of the statistical procedure for grid refinement, many 
elements are stuffed in the large error zone and this leads to a 
rapid decrease in the values of the maximum as well as the mean 
errors. Indeed, since the statistical procedure is applied on 






0.0025 
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successive meshes several times until the error tolerance criteria 
are met, within one time step the maximum and mean errors are 
reduced by a few orders of magnitude. Moreover, a natural 
gradation of grid spacing results which is well in accordance with 
the solution activity and consequently the local residue 
distribution approaches zero. In the moving front region, the 
mean error undergoes mild fluctuations due to the fact that there 
is some time lag between front movement and grid adjustment; in 
other words, the grid-adaption occurs only after sufficient 
residue has accumulated. This inference is supported by the trend 
that the maximum error increases with time sharply, reaches a peak 
value and then decreases rapidly due to grid-adaption. On the 
converged mesh, the maximum error is close to the mean, indicating 
an effective equi- distribution of error by the statistical 
technique (Fig. 4.13a). 

(d). Adaptive mesh : 

Figs . 4 . 14a and 4 . 14b show the adapted meshes at various time 
instants in the h and hp refinement procedures. The final 
adaptive meshes have been plotted after the convergence criterion 
is satisfied. Also, the number of elements have been indicated 
in the side tables presented along with the figures. The 
computations have been started with a uniform mesh of N = 10 
elements. Further, the results of the h and h-p algorithms have 
also been presented in the tables (4.4a,b) at different time 
points. 

The figures clearly illustrate the ability of the adaptive 
scheme to concentrate grid points around the characteristic wave 
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Time 

Adp 

Elements 


Nodes 

Apos .err 

clni . err 

Eff.inc 

Rem 

Ref 

Stf 

0.00005 

1 

0 

0 

0 

10 

11 

2.9327 

2.0 () () 3 

1-4193 


2 

2 

1 

4 

11 

12 

0.6807 

1 .8(54 9 

0.3650 


3 

0 

1 

4 

14 

15 

0.3248 

1 . 8 500 

0 . 1756 


4 

3 

3 

8 

16 

17 

0. 1357 

1 . 8 2 6 6 

0.0743 


5 

3 

2 

8 

19 

20 

0.0557 

1 . 8 2 9 7 

0.0305 

0.00010 

1 

0 

0 

0 

19 

20 

0.1374 

1 . 6 7 2 1 

0.0822 


2 

0 

4 

10 

25 

2 6 

0.0708 

1 . 6 7(H) 

0.0424 

0.00045 

i 

0 

0 

0 

34 

35 i 

0.0828 

1 .0297 

0 . 0804 


2 

0 

9 

22 

47 

48 

0.0422 

1.0288 

0.0411 

0.00110 
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0 

0 

0 

47 

4 8 
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5 

8 

18 

52 

53 

0 . 04 7 6 
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0.08 3 3 
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0 
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52 

53 
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24 

64 
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134 

0 . 0762 

0. 109 4 

0.6960 


2 

22 

19 
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0 
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0 
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1 

0 

0 

0 
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154 
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2 

39 

24 

58 
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0.0373 

0.07 7 5 

0.4832 
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1 

0 

0 

0 
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1 

0 

0 

0 
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0.0725 
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1 

0 

0 

0 
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0.0364 
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0.6204 

0.07880 

1 

0 

0 

0 
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13 3 

0.0613 

0.0744 

0.8240 

0.08380 

1 

0 

0 

0 

140 
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0.0761 

0 . 0847 

0.8984 


2 

29 

24 

60 

147 
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: 0.0342 

0.0769 

0.4447 

0.08405 

1 

0 

0 

0 

147 
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0.0347 

0.0701 

0.4942 

0.09080 

1 

0 

0 

0 

153 
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0.0766 

0.0867 

0.8833 


2 

36 

24 

58 
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0.0333 

0.0784 

0.4245 

0.09105 

1 

0 

0 

0 
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I 0.0349 

0.0710 

0.4916 

0.09305 

1 

0 

0 

0 
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2 

36 

24 

58 

149 
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0.0322 

0.0783 

0.4110 

0.09880 

1 

0 

0 

0 
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0.0298 

0.0496 

0.6010 

0.09955 

1 

0 

0 

0 
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0.0248 

0.0440 

0.5638 


Table 4.4a. Test problem results for h-version strategy 



163 


Time 

Adp 

Deg 

Rem 

N 

Dof 

IEHSB5I 

BnVtB 

Ef f . ind 

0.00005 

1 

1 

0 

10 

11 

2.9327 

2.0663 

1.4193 


2 

1 

2 

8 

9 

2.9327 

2.0663 

1.4193 


3 

2 

0 

8 

10 

1.3210 

2.0599 

0.6413 


4 

3 

0 

8 

12 
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5 
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Table 4.4b. Test problem results for h-p version strategy 
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x - Ct - 0, where all the activity takes place. In the h 
refinement procedure , the adaptive algorithm automatically 
removes , retains and stuffs nodes simultaneously. It is however , 
possible that some unnecessary nodes exist in regions where no 
solution activity is taking place. In the present work, removal 
of nodes has been carried out if the error of the neighbouring 
elements goes below a minimum error tolerance of 10~ 6 . Some times 
this tolerance value is very much below the mean error and so even 
for regions with no solution activity, node removal does not take 
place. In order to eliminate this problem, dynamic setting of the 
minimum tolerance limit ( as a fraction of the mean error ) can be 
done, as discussed in Chapter 5. On the other hand, over 
refinement of the grid occurs some times due to the statistical 
refinement procedure. In this scheme, if the local error is 
greater than the mean, element refinement is performed. At later 
stages of grid adaption, the mean error is very small and it is 
prone to be affected by round off error. Thus, at few instants, 
elements which need not be refined are also divided further. 
Occasionally, excessive refinement may lead to deterioration of 
the computed solution due to domination of round off errors; but 
this can be easily checked by specifying a minimum grid spacing. 
These difficulties are, however, rare exceptions and as seen from 
the Figs. 4.14a-b, the performance of the h and h-p refinement 
strategies is remarkably good. The higher degree approximations 
of the p version tend to produce oscillations in the solution. 
Therefore, the local elemental errors may not behave well, 
especially adjacent to the moving front. Due to these 
oscillations, some nodes are not properly removed from low 
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gradient regions during hp adaption- From both Figs. 4.l4a-b, it 
is clear that the nodes are distributed in anticipation of the 

front motion. 

Comparing the relative performance between the h and hp 
refinement methods, the following observations can be made. i n 
general, the h-version is more efficient in bringing down the 
error than the hp- version. Therefore, both the maximum and mean 
error values are lower in the h-version (Fig. 4.13a and 4.13b). 
Further, after each grid adaption, the maximum error is brought 
very close to the mean error by the h refinement strategy, 
implying a very effective equi-distribut ion of error. On a 
comparative scale, the degree of equi-distribut ion by the hp 
refinement method is not very good. Also, the error decreases 
slowly when the degree of interpolation is increased by the 
statistical procedure. In the present work, if the maximum 
allowed polynomial degree is reached without reaching error 
convergence, the hp -version mesh is converted into an 
approximately equivalent h-version mesh. After this conversion, 
the error may slightly increase or decrease. Due to the slow 
convergence, the hp method requires more number of adaptions. The 
higher values of the maximum and the mean errors observed in the 
hp version can be attributed to the following reasons. It is 
observed that for interpolation order greater than 3 , the computed 
solution is more prone to numerical noise, which is seen in the 
form of oscillations in the computed solution. Also, as the 
degree of interpolation is increased in one element, the computed 
error of the neighbouring elements tends to increase. On various 
accounts, the performance of the h-version refinement appears to 
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be superior , although the number of elements required by the hp 
version is less in order to meet the same error criterion. The 
computational effort per adaption is also less in the h-p 
version due to the application of the hierarchical basis 
function concept. 

(e) . Effect of solution smoothening on true gradient error 

estimation : 

The effects of the smoothening procedure upon the error 
measures are shown in the Figs. 4.15a-c. Due to the linear degree 
approximation used in the h -version procedure, the computed 
gradient is constant within each element. On the other hand, the 
exact gradient varies from point to point even within an element. 
For this reason, piece-wise cubic spline smoothening procedure has 
been applied to the computed solution in each element, for 
improving the estimates of true gradient error and the effectivity 
index. For smoothening purposes, the average value between the 
forward and backward derivatives of the gradient at the nodes are 
first evaluated. Using these average gradients as well as the 
values of the computed solutions at the nodes, piece-wise cubic 
Hermite interpolating polynomials have been employed for 
smoothening. 

It is seen that the true gradient and relative errors are 
somewhat reduced after smoothening, thereby bringing the value of 
the effectivity index closer to unity for the non-uniform mesh. 
However, the magnitudes of the oscillations in the true gradient 
error and the effectivity index increase after smoothening which 
is not an acceptable feature. In any case, the previously 
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observed trend that the a posteriori error is slightly less than 
the true gradient error for non— uniform meshes , holds with 
smoothening also. For the hp refinement method, due to its 
inherently smooth nature of solution (of higher degree 
approximating polynomials) , the smoothening of the computed 
solutions has not been performed while estimating true gradient 
error for the test problem. 

(III). Grid speed : 

On a coarse mesh the grid speed calculation is inaccurate. 
So, we consider the solution at a particular time when adequate h- 
refinement has been achieved and a fully established front moves 
at a uniform velocity. 

Fig. 4.16. shows the behaviour of solution profile at various 
times during the uniform front propagation stage. The exact 
solution and the solution by grid speed have been superimposed. 
It is observed from the figure that the computed solution has good 
correspondence with the exact solution. Because of the highly 
non- uniform mesh which has been considered here, the solution 
front, moves in a smooth way. However, some difficulties do arise 
while predicting the solution in a small neighbourhood near the 
boundaries x = 0.0 and x = 1.0. The reason for this is that if 
all the grid points move with same velocity, collapsing and 
stretching of elements occur near the boundary which affect the 
solutions. In order to avoid this, special treatment of grid 
movement is necessary near the domain boundaries. 

In the computational experiments, the boundary nodes of the 
domain have been locked and the interior nodes are allowed to move 



according to uniform grid speed. After a few t i mo iterations, the 
nodes which are near to the boundary may shoot outside the 
solution domain. These nodes which cross t ho boundary are 
reintroduced at the beginning of the mesh, as soon in the 
Fig. 4.17. Although, it is not necessary to reintroduce these 
nodes in the present problem, for more general cases it may be 
necessary. It is clear from the Fig. 4.17. that the nodes are 
distributed according to the front. The computed grid velocity 
for this problem matches very well with the exact 1 ront speed 
value of 10. 

Fig. 4.18a. shows the variation of a posteriori error with 
time in the wave front region. It is observed that the a 
posteriori error increases very slowly, but lor the entire wave 
propagation time, the error is very small (smaller than the 
tolerance used for h and hp methods). This shows that the 
elemental error indicators also move according to 1 ront . The 
small magnitude of the errors also implies that the chosen initial 
mesh( obtained by h- version ) is adequate for resolving the sharp 
spatial gradients. The slow increase in error with respect to 
time may be attributed to a small lag between the grid movement 
and the front movement. 

A similar behaviour is observed in the Fig. 4 . 18 b in which 
the true gradient error has been plotted. The decrease in the true 
gradient error at initial time shows that the gradient of the 
space discretization error decreases rapidly due to very fine mesh 
around the front location at that time. At later times, since the 
grid movement does not exactly coincide with the front, true 
gradient error increases slowly. The relative error behaviour 
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shown in Fig. 4.18c is same as that of the gradient error. 

Results for the effectivity index versus time are shown in 
Fig. 4 . 18d . From the graph, it can be seen that the effectivity 
index increases rapidly at initial times and it remains almost 
constant till the final time is reached. Its value is not too far 
from unity, indicating a proper relationship between the a 
posteriori and true gradient error estimates. From Figs. 4.18a-d, 
it is clear that the initial variations of these error measures 
for grid-speed method depend on the starting solution and the 
mesh. For later times, the variations are governed by the 
difference between the grid and the front speeds. If the starting 
solution is good and the estimate of the grid speed is accurate, 
error measures can be maintained below the tolerance level for the 
whole duration of the problem. It is, however, advisable to check 
whether the solution obtained from grid speed procedure satisfies 
error criteria or not, from time to time. Thus, the grid speed 
technique can be an inexpensive way of constructing adaptive 
meshes for problems with constant front speeds. 

4.8.2. Application of adaptive strategy to the flame propagation 
problem : 

In this section, results obtained by h-version, hp-version 
and grid-speed methods for the flame propagation problem described 
by equations (4.4) are presented. The solution profiles 
(temperature and reactant density) and the variation of error 
measures have been plotted at various time instants, for both 
uniform and non-uniform (adaptive ) meshes. The evolution of 
adaptive meshes with time have also been shown for h, h-p and the 
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grid speed strategies. 

(I). Uniform mesh results : 

For the non -linear flame problem considered in the present 
study, exact solutions are not available. There lore, it is not 
possible to estimate the true gradient error in the computed 
solution. However, in the previous section, the behaviour of the 
true gradient and the a posteriori errors wore shown to be 
identical for the test problem which has similar characteristics 
as the flame problem. Therefore, a posteriori error which is 
computable from the elemental residues, is taken as the error 
measure for grid- adaptions. Bieterman and Babuska ( 1982a-b; 

1986) have also used the same error estimate for the problems of 
this nature. The flame problem was solved on a sequence of 

uniform meshes with 40, 80, 160 and 320 elements. The reasons for 
considering the uniform meshes have been cited in the previous 
section. For the flame problem also, uniform grid results can be 
used as reference, since they are free of the moving grid effects. 

(a). Temperature and density profiles : 

In Figs. 4 . 19a-b and 4 . 20a-b, the temperature and density 
profiles for the flame propagation problem have been shown for 
uniform meshes with 40, 80, 160 and 320 elements. For the coarse 
mesh ( N = 40 ) , the maximum flame temperature is unrealistically 
high; besides, oscillations are observed at the edges of the flame 
front due to inaccurate resolution of the gradients. In fact, as 
described in detail by Dwyer (1983b), the flame profile reaches 
its equilibrium shape some time after ignition, due to the balance 



! 1 
\ ] 

!! 

/ i 


it 


i 1 
} , 

! i 


UL 


12 


0.4 


0.6 



_J , I I I L 

0.2 0.4 0.6 

X - SPATIAL DISI 

FLAME PROBLEM (Uniform Mes 




Mesh 

, ! 
Line type 

N = 40 

N = 80 

— 


Act’ 

Time 

Eles I 

Nods 

1 

0 

.000010 

40 

41 

2 

0 

.000150 

40 

41 

3 

0 

.000300 ! 

40 

41 

4 

0 

. 000800 : 

40 

41 

5 

0 

.001550 

40 

41 

a 

0 

.002550 

40 

41 

7 

0 

.003550 

40 

41 

8 

0 

.004550 

40 

41 

9 

0 

.005050 

40 

41 

10 

0 

.006050 

40 

41 


Adp 

Time 

Eles 

Nods 

1 

0 

.000010 

80 

81 

2 

0 

.000150 

80 

81 

3 

0 

.000300 

80 

81 

4 

0 

.000800 

80 

81 

5 

0 

.001550 

80 

81 

6 

0 

.002550 

80 

81 

7 

0 

.003550 

80 

81 

8 

0 

.004550 

80 

81 

9 

0 

.005050 

80 

81 

i 

10 

0 

.006050 

80 

SI 


) : SOLUTION PROFILE 








temperature 



Fig. 4.20 


1 

i\ t 

v\ \\ 


0.2 0.4 0.6 0.8 1.0 1.2 

X - SPATIAL DISTANCE 

FLAME PROBLEM (Uniform Mesh ) : SOLUTION PROFILE 




Mesh 

Line type 

N = 160 

N = 320 

— 


Adp 

Time 

Eles 

Nods 

1 

0 .000010 

160 

161 

2 

0 .000150 

160 

161 

3 

0 .000300 

160 

161 

4 

0 .000800 

160 

161 

5 

0 .001550 

160 

161 

6 

0 .002550 

160 

161 

7 

0 .003550 

160 

161 

8 

0 .004550 

160 

161 

9 

0 .005050 

160 

161 

10 

0 .006050 

160 

161 


Adp 

Time 

Eles 

Nods 

1 

0 .000010 

320 

321 

2 

0 000150 

320 

321 

3 

0 .000300 

320 

321 

4 

0 .000800 

320 

321 

5 

0 .001550 

320 

321 

6 

0 .002550 

320 

321 

7 

0 .003550 

320 

321 

8 

0 .004550 

320 

321 

9 

0 .005050 

320 

321 

10 

0 .006050 

320 

321 








177 


between the reaction and the heat conduction process. Since, the 
heat conduction process depends on the solution gradient, it is 
understandable that an inaccurate estimation of the 
temperature/density gradients will lead to an inaccurate 
prediction of the reaction rate and the flame temperature. The 
overshoot in flame temperature is eliminated for finer meshes and 
the oscillations in the solution profiles are also absent which 
can be seen from the Figs. 4.2 0a and 4.20b. For the coarse mesh, 
although the solutions are highly inaccurate, all the flame 
features such as equilibrium flame shape, constant speed of front 
movement etc. , are observed. For more refined meshes (with 
N = 80, 1 60 or 320), converged equilibrium flame shapes are 
obtained. Just after ignition, the flame profile exhibits a small 
hump in its shape. The presence of the hump is attributed to the 
shift in equilibrium shape due to the higher heat conduction 
effect in the proximity of the wall (x = 1.0 ). When the flame 
moves away from the wall, the hump in the flame profile 
disappears . 

Regarding the density profiles, the following observations 
are made. For coarse mesh ( N = 40 ) , the density profile also 
has fluctuations in the form of overshoots or undershoots near the 
edges of the flame. Indeed, after the flame front has moved 
through a particular location, the density of the reactants at 
that location is expected to go zero due to complete burning. 
However, for the coarse mesh, there is a tendency to predict small 
negative values of density in the burn-out region, due to 
numerical oscillations. In the computational scheme, these non 
physical negative values have been forcefully made equal to zero. 
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The reaction term acts as a heat source in the temperature 
equation, while it is a sink term for the reactant density 
equation. For this reason, the density profile has just the 
inverse shape of the temperature profile. All the other features 
such as front movement etc., are depicted by the density profile 
also. For the number of elements N = 160 and 120, identical 
profiles are obtained, for each instant of time. 

In general, for N = 40, the convergence of the temperature or 
density profiles posed some difficulties during iterative method 
of solution due to large fluctuations in the predicted values. For 
fine meshes, this convergence problem did not arise. The computed 
shapes and front locations widely differ between the meshes of 
N = 40 and 80 while excellent agreement is observed between the 
predictions of N == 160 and 320. It i.s noted that, the front speed 
predicted by the coarse mesh is also highly inaccurate. 

(b). The behaviour of a posteriori error estimate : 

In Figs. 4.21a-d, the variation of the a poster ior i error and 
the solution gradient with time are shown tor different number of 
elements. The trends of a posteriori error seen for the flame 
problem are very similar to those of the test problem. During 
the ignition/ extinction phases of the flame, the a posteriori 
error is small due to the low value of solution gradients. An 
interesting feature seen for the coarse ( N = 4 0 ) mesh is that 
the flame extinguishes earlier compared to the finer meshes, due 
to higher flame velocity. In fact, computations have not been 
carried till extinction for N = 80, 160 and 320. 

In the flame propagation region systematic high amplitude 
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fluctuations are observed in the error, similar to the results of 
the test problem with uniform elements. For N = 40, the a 
posteriori error exhibits large oscillations and the amplitude of 
the fluctuations are very high. A similar trend has been observed 
for N = 80, although the magnitude of the fluctuations have been 
reduced . 

As the number of elements increases, the amplitude of the 
fluctuations decreases and for N > 160, it is negligibly small. 
It is observed that for each mesh the maximum error occurs 
slightly after ignition when the hump formation is seen in the 
solution profile. This is due to the fact that the nonlinear 
combustion reaction is most active at that particular time, 
possibly because of high temperature occurring in the reactant 
mixture. When the flame front moves with a constant velocity, the 
process of combustion reaction, heat transfer and front movement 
are in perfect equilibrium with each other; thus the mean value of 
the a posteriori error reaches a steady value in this region. 
However, the magnitude of the this steady value of error, 
decreases with mesh size. Analogous to the results of the test 
problem, the error is seen to be linearly proportional to the mesh 
size, especially, for the fine meshes. For coarse meshes, the 
error seems to decrease faster than h, when the grid is refined. 
It is to be noted from the Figs. 4.21a-d that even on the fine 
mesh ( N = 320 ) the a posteriori error does not meet the error 
tolerance criterion. 

Thus, from the above discussions, it is clear that the error 
behaviour for the flame propagation problem has very similar 
characteristics as that of the test problem. Whatever differences 
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exist between the trends of these two problems, can be attributed 
to the inherent nonlinearity of the flame problem and also the 
differences in the front width as well as the front speed between 
the two problems. 

(II). Non-uniform mesh results : 

The analysis carried out by the uniform grids shows that a 
large number of grid points are necessary to resolve the flame 
front which satisfies the error tolerance criterion. So, the 
solution procedure becomes computationally very expensive when 
uniform grids are used. In order to have an accurate physical 
simulation, adaptive grids are necessary which automatically 
identify the location of the active zone. Results obtained with 
h, hp and grid speed strategies are presented in the ensuing 
sections. 

(a). The temperature and the density profiles : 

The solution profiles obtained by h version and the uniform 
mesh N = 32 0, have been superimposed in the Figs. 4.22a-b. 
Further, a comparison between the solutions obtained by the h and 
hp version has also been shown the Figs. 4.23a-b. In both the 
figures, the solutions were obtained after satisfying the same 
tolerance limit on the a posteriori error. The h and hp 
computational procedures were started with a coarse mesh of 10 
uniform elements. The profiles predicted by the h version and 
the N = 320 meshes are similar. However, there is a small 
difference in the locations of the flame fronts indicating a 
slightly inaccurate prediction of flame speed by the uniform 
mesh. Such a difference is expected since the uniform mesh does 
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not employ the mesh points optimally and also it door, not satisfy 
the error criterion. The solutions predicted by the h and the h-p 
strategies are identical at all time instants, thereby verifying 
the accuracy of each other. The hp method requires less number of 
elements for the same computable accuracy, but it involves more 
adaptive cycles. The above comparisons establish the .‘supremacy of 
the error based adaptive solution strategies of over the fixed 
grid numerical computations. The important point to be noted here 
is that for the uniform element solution one needs to know the 
number of elements which will produce solution profiles up to the 
desired accuracy; for the h and hp methods, the number of elements 
is automatically decided - but the adequate tolerance, level is 
to be found by trial and error. In the error based methods, 
however, one has direct control over the solution error and the 
grid distribution is also more rational. 

(b). The behaviour of a posteriori error : 

In Fig. 4.24, the variations of a posteriori error -and the 
norm for the computed gradient h and hp grid ref inement, strategies 
have been plotted. The initial number of elements have been taken 
as 10 in both cases. Because of the coarse initial mesh, the 
error shoots up to a very high value in the first time step for 
both methods; but due to the error based adaption strategy, it is 
brought down rapidly as well, within a single time stop. As 
ignition proceeds, the a posteriori error slightly increases and 
reaches a peak value with the occurrence of the hump in the 
temperature profile. As explained in the case of uniform meshes, 
this peak in the error is linked to the highest rate of reaction 
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prevailing at this juncture. As soon as the t lame attains the 
equilibrium phase of uniform propagation, the a posteriori error 
fluctuates, with a steady mean value. The decrease in the a 
posteriori error seen in the figure after each adaption at any 
time is due to the effect of the statistical procedure which 
reduces the local residue in the greater solution activity region. 
Also, the L 2 norm of the computed gradient is constant in the wave 
front region except at initial stages for both the methods . 

Even though the necessary mathematical theory har. not been 
developed for the implementation of the hp method for the flame 
propagation problem, the computational experience shows that the 
present a posteriori error estimate for the linear problems works 
well. The number of adaptions required for the hp strategy is 
larger due to slower sensitivity of error to grid adaption, as 
compared to the h version procedure . Hut with such minor 
differences excluded, both the methods perform very wolf in 
controlling the solution error within the spoci f iod • t olerance. 
Most of the trends exhibited by the a posteriori error during the 
application of the adaptive h and hp strategies, tor the flame 
propagation problem bear similarities with the uniform mesh 
results of the flame propagation as well as the test problem 
results with uniform and non-uniform grids. This indicates that a 
posteriori error which is based on the residue, is a rugged error 
measure suitable for grid adaption purposes in both 1 inoar and 
non-linear problems involving front propagation. Home of the 
numerical results obtained by the h and hp refinement methods for 
various time instants have also been presented in the tables (4.5) 
and (4 . 6a-b) . 
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Time 


Elments 





Adp 

Rem 

Ref 

Stf 

N 

Nodes 

Apos . err 

Com. grd 

0.000010 

1 

0 

0 

0 

10 

11 

24.52254 

0.2078 


2 

1 

1 

4 

12 

13 

3 . 06718 

0.4067 


3 

2 

1 

4 

13 

14 

0.40373 

0.6596 


4 

5 

1 

4 

11 

12 

0.08006 

0.7418 


5 

3 

2 

6 

12 

13 

0.03707 

0.7495 

0.000020 

1 

0 

0 

0 

12 

13 

0.11114 

1.1890 


2 

0 

2 

6 

16 

17 

0.11693 

1.1835 


3 

0 

3 

10 

23 

24 

0.05885 

1.1824 

0.000060 

1 

0 

0 

0 

31 

32 

0.07549 

2.6877 


2 

2 

8 

20 

41 

42 

0.06521 

2.6872 

0.000115 

1 

0 

0 

0 

59 

60 

0.07878 

4.4325 


2 

3 

15 

36 

77 

78 

0.05337 

4.4327 

0.000145 

1 

0 

0 

0 

113 

114 

0.07601 

6.4919 


2 

1 

23 

52 

141 

142 

0.05562 

6.4918 

0.000160 

1 

0 

0 

0 

141 

142 

0.08230 

8.1623 


2 

1 

26 

58 

172 

173 

0.06618 

8 . 1623 

0.000200 

1 

0 

0 

0 

232 

233 

0.07582 

10.0157 


2 

32 

72 

186 

314 

315 

0.03995 

10.0156 

0.000245 

1 

0 

0 

0 

314 

315 

0.07713 

9.8579 


2 

72 

73 

182 

351 

352 

0.03751 

9.8578 

0.000360 

1 

0 

0 

0 

352 

353 

0.07791 

9.2782 


2 

95 

54 

144 

347 

348 

0.03911 

9.2781 

0.000500 

1 

0 

0 

0 

346 

347 

0.07586 

8.9779 


2 

95 

77 

196 

370 

371 

0.03682 

8.9777 

0.000855 

1 

0 

0 

0 

349 

350 

0.07583 

8.7898 


2 

93 

57 

138 

337 

338 

0.04011 

8.7896 

0.001165 

1 

0 

0 

0 

360 

361 

0.07806 

8.7668 


2 

97 

57 

150 

356 

357 

0.03791 

8.7666 

0.001670 

1 

0 

0 

0 

351 

352 

0.07850 

8.7621 


2 

95 

71 

146 

331 

332 

0.04474 

8.7619 

0.002635 

1 

0 

0 

0 

367 

368 

0.07552 

8.7614 


2 

99 

60 

146 

354 

355 

0 . 03883 

8.7612 

0. 003355 

1 

0 

0 

0 

333 

334 

0.07675 

8.7616 


2 

90 

75 

154 

322 

323 

0.04212 

8.7614 

0.004855 

1 

0 

0 

0 

380 

381 

0.07877 

8.7623 

; 

2 

108 

54 

122 

340 

341 

0.04269 

8.7621 

0.005810 

1 

0 

0 

0 

347 

348 

0.07504 

8.7620 


2 

93 

65 

164 

353 

354 

0.03976 

8.7618 

0.005960 

1 

0 

0 

0 

351 

352 

0.07687 

8.7621 


2 

97 

62 

156 

348 

349 

0.03796 

8 .7619 

0.006035 

1 

0 

0 

0 

348 

349 

0.07728 

8.7625 


2 

96 

65 

168 

355 

356 

0.03686 

8 .7623 


Table. 4.5. h version results for flame propagation problem 
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Time 

Adp 

Deg 

Rem 

N 

Dof 

Apor. .err 

Com.grd 

0.000010 

1 

1 

0 

10 

11 

24 . 52264 

0.2078 


2 

1 

1 

9 

10 

24 . 52254 

0.2078 


3 

2 

0 

9 

11 

8.1734 8 

0.2881 


4 

3 

0 

9 

12 

3 . 964 37 

0 . 3722 


5 

4 

0 

9 

13 

2.23357 

0.4538 


6 

5 

0 

9 

14 

1.34661 

0 . 5292 


7 

1 

0 

12 

13 

0.48328 

0 . 6393 


8 

1 

3 

9 

10 

0.48328 

0 . 6393 


9 

1 

2 

7 

8 

0.48328 

0 . 6393 


10 

2 

0 

7 

9 

0. 1 5986 

0.7236 


11 

3 

0 

7 

10 

0.09315 

0.7465 


12 

4 

0 

7 

11 

0 . 08333 

0.74 91 


13 

5 

0 

7 

12 

0.07093 

0.7493 

0.000040 

1 

1 

0 

25 

26 

0.08108 

1 . 9831 


2 

1 

2 

23 

24 

0 . 08 2 06 

1.9832 


3 

2 

0 

23 

30 

0.05484 

1 . 9834 

0.000110 

1 

1 

0 

43 

44 

0.08978 

4.2384 


2 

2 

0 

43 

56 

0 .1 3 3 4 1 

4 .2382 


3 

3 

0 

43 

57 

0 . 1206 5 

4.2382 


4 

4 

0 

43 

58 

0.1132 6 

4 .2382 


5 

5 

0 

43 

59 

0. 1 0736 

4 . 2 382 


6 

1 

0 

47 

48 

0 . 08268 

4 .2383 


7 

1 

1 

46 

47 

0 . 08269 

4 . 2383 


8 

2 

0 

46 

60 

0 . 06729 

4 . 2382 

0.000255 

1 

1 

0 

236 

2 37 

0 . 07992 

9 . 7 9 3 0 


2 

1 

10 

226 

227 

0 . 080 1 1 

9 . 7 9 3 0 


3 

1 

1 

225 

22 6 

0 . 08 0 1 5 

9 . 7930 


4 

2 

0 

225 

278 

0 . 0 5 1,3 9 

9 . 7930 

0.000790 

1 

1 

0 

211 

212 

0.079 1 8 

8 . 8032 


2 

1 

10 

201 

202 

0.0794 1 

8 . 8032 


3 

1 

1 

200 

201 

0.07945 • 

8 . 8032 


4 

2 

0 

200 

246 

0.05102 

8 . 8032 

0.001270 

1 

1 

0 

221 

222 

0.07489 

8 . 7 64 8 


2 

1 

14 

207 

2 08 

0.07628 

8.7647 


3 

1 

1 

206 

207 

0.07531 

8.7647 


4 

2 

0 

206 

268 

0.04535 

8.7648 

0.001500 

1 

1 

0 

220 

221 

0.07531 

8.7622 


2 

1 

7 

213 

214 

0.07545 

8.7622 


3 

1 

1 

212 

213 

0.07547 

8.7622 


4 

2 

0 

212 

263 

0.04794 

8 . 7623 

0.001720 

1 

1 

0 

217 

218 

0.07503 

8.7612 


2 

1 

12 

205 

206 

0,07529 

8.7612 

0.002020 

3 

2 

0 

205 

256 

0.04864 

8.7612 

1 

1 

0 

209 

210 

0.08069 

8.7614 


2 

1 

11 

198 

199 

0.08095 

8.7614 


3 

2 

0 

19 8 

249 

0.05474 

8.7614 


Table. 4.6a. h p version results for flame propagation problem 
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(c) . The effect on Maximum and Mean of local error indicators : 

The plots of maximum and mean local error variations with 
respect to time are presented in Figs. 4.25 a-b, for the h and hp 
methods. The difference between the maximum and the mean errors 
quantifies the extent of adaption required at the location of the 
maximum error. At initial times ( during the ignition phase), the 
difference between the maximum and the mean errors is high, 
indicating the great need for grid adaption. Indeed, at this 

time, the mesh structure changes from coarse to line by rapid 
stuffing and removal of elements. After the completion of the 
ignition phase ( and the hump region ), the difference between the 
maximum and the mean errors reduces considerably, indicating a 
well graded mesh. Occasionally, due to slight mismatch between 
the existing mesh and the front location, the maximum error 
exhibits high peaks. The mean local error, on the other hand, has 
a steady value during equilibrium propagation phase. The peaks in 
the maximum error are followed by sharp falls due to the grid 
adaption process which sets in at these peak points. In fact, 
after each adaption the maximum error is brought close to the mean 
value by equidistributing the error, which indicates the efficient 
application of the statistical procedure. Due to the high 
sensitivity of the h version, the difference between the maximum 
and the mean becomes very small after only one grid adaption. 
For the hp method, in fact, the statistical procedure is not as 
effective as in the h method. The number of adaptions required in 
hp version required is also much more as discussed earlier. Due 
to this large number of grid adaptions and also the conversions 
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from p to h-version meshes, the mean for the hp error forms a 
thick band with respect to time. In general, both the mean and 
the maximum of local errors for the h method are below those of 
the hp method. Most of these features of local errors for the 
flame problem are similar to the trends observed lor the test 
problem when the non-uniform meshes are used. 

(d). Behaviour of local error indicators and distribution : 

The distributions of the square of local error indicators at 
various time instants and adaption levels have been plotted in 
Figs. 4.26a-f for the h refinement method. The time levels have 
been chosen such that the ignition, flame propagation and 
extinction phases of the flame are all covered. The computations 
have not been performed till extinction is complete. With a 
starting mesh of 10 uniform elements, the error indicator is very 
large in magnitude adjacent to the boundary where the flame 
ignites ( see Fig. 4.26a ). As the mesh is ret ined according to 
the statistical procedure in the large error zone, the error 
indicator decreases rapidly. Also, in the hump region the of 
solution (i.e. near x = 1.0 boundary at t - 0.0003), the 
behaviour of the local error indicators exhibit fluctuations as 
seen from Fig. 4.26b. By examining the ratios of the maximum 
error and the corresponding step length the decrease in error 
after grid adaption is observed to be proportional to h. It is to 
be noted, however, that the precise dependence of the error 
indicators on the local step length is very difficult to ascertain 
due to the fact that the predominant contribution to the local 
residue arises from the non-linear source term. This aspect is 
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addressed in greater detail in Chapter 5. in Fig. 4.26a, it is 
seen that the change in the number of elements after each adaption 
is very small, although the local error has decreased by orders of 
magnitude. The reason for this is the simultaneous removal of 
nodes from unwanted zones, along with grid refinement. 

At different time instants, the local error distribution 
exhibits similar features. The maximum value of the local error 
occurs at the location of maximum reaction rate ( flame front 
region) . Occasionally, due to excessive coarseness of the mesh, 
the local error increases for some adaptions in the non-reaction 
zones also. However, this increase is confined to one or two 
elements only and gets immediately rectified by proper 
grid-spacing adjustment in the subsequent adaptions ( see. Fig. 
4.2 6b-c ). Thus, with regard to the satisfaction of the global 
error criterion, these large elements do not have significant 
influence. In general, the magnitude of the maximum local error 
is brought down considerably ( by almost an order of magnitude ) 
by each h version grid-adaption. 

The width of the flame zone can be correlated with the region 
where' close peaks in the local error are observed. It is clear 
from the Figs. ( 4.2 6 a-f ) that during the ignition phase the 
flame region is very thin but its width grows with time and 
reaches a constant value during the equilibrium propagation phase. 
Due to this increase in width, more number of elements contribute 
to global error after the ignition phase. Therefore, if the same 
global error value is maintained for all times, the satisfaction 
of error tolerance criterion leads to very small values of the 
local error for later time instants, than in the beginning. This 
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effect is clearly observed in the figures. For time values 
greater than 0.0003, the fully developed equilibrium flame front 
is established and the number of elements does not vary much 
during this period with respect to time as well as adaptions 
within a time step. The main change that occurs is the 

redistribution of dense grid point zone so as to match the 
shifting location of the flame front region. Because of this, the 
band of local error peaks observed in the Figs. 4.26a-f, shifts 
its location with respect to time level or adaption within a time 
step. As regards the number of adaptions required for satisfying 
the error criterion, it is high only for the first few time steps 
and for later times, the error tolerance is met within 2 
adaptions, since the grid distribution becomes almost optimal at 
every stage and deviations from optimality are minimal. 

The mesh distribution at different time instants and adaption 
steps, are shown in Figs. 4.27a-b for the h version method. The 
solution procedure was started with an initial mesh of 10 uniform 
elements. In the first time step, simultaneous refinement and 
element removal (derefinement) are seen during each adaption. It 
is observed that the refinement procedure is somewhat slower than 
the refinement. There are two reasons for this :(1) In the 

present algorithm when two neighbouring elements have error below 
a prescribed minimum tolerance level, they are collapsed into a 
single element. In a given adaption step, this procedure is 
carried out on every pair of existing elements satisfying the 
above criterion. However, for the new elements created by the 
merger of two elements in the existing mesh, the check for element 
collapse will be performed only in the next adaption. (2) The 
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magnitudes of the error indicators and their statistical measures 
( mean, standard deviation, etc., ) decrease with grid adaption. 
Thus, elements which would not have qualified for collapse become 
so in the subsequent adaption steps. 

It should be mentioned in passing , however, that the adaption 
process on the coarse mesh is influenced by refinement or 
derefinement occurring anywhere in the domain; therefore, 
execution of removal mechanism for elements in the non reactive 
region depends on the grid adaption occurring in the other places. 
These features are clearly seen in the evolution of the adaptive 
mesh at t = 0.00002. With respect to time, it is observed that the 
numerical grid faithfully follows the front movement. In fact, if 
the element distribution is not optimal, then within the same time 
step the mesh is adjusted appropriately in subsequent iterations. 
It is observed in the figures that during this iteration process, 
the dense mesh region shifts to the exact location of the flame 
front. The width of the dense mesh region which is a measure of 
the thickness of the flame front, remains more or less constant 
during the equilibrium flame propagation-phase. It is observed 
that the total number of elements in the solution domain 
fluctuates with time, even in the equilibrium phase. The reason 
for this is that the tolerance criterion for error is an 
inequality constraint only and therefore, any number of elements 
above a certain minimum value will satisfy the error tolerance. 

In Figs. 4.28a-f, the distribution of the squares of the 
local error indicators have been shown for the hp refinement 
procedure. The corresponding adaptive meshes and the degree of 
interpolation used within the elements are depicted in 
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Figs. 4.29 a-c. The number of elements and the total degree of 
freedom in the problem have also been indicated for each instant 
of time and adaption level. Here, the grid adaption process has 
been divided into three steps, namely : (i) Derefinement where 
local error is below the minimum tolerance, (ii) Increase of 
interpolation degree in elements where the local error is above 
the mean value, (iii) conversion of p - version mesh into an 
approximately equivalent h mesh, for all the elements which have 
attained the maximum degree of interpolation permissible without 
meeting the tolerance criterion. For each time step, this 
procedure has been performed sequentially, until the global error 
criterion is met. During computation, each of these procedures 
have been carried out in a separate adaption step, for the sake of 
convenience. However, in the figures, only a few representative 
adaption steps have been shown for each time step and the numbers 
indicated do not correspond to the exact adaption step number. 

It is observed from Figs. 4.28a and 4.29a that the local 
error decreases rather slowly when the degree of approximation is 
increased. Further, the decrease in the local error in some zones 
may lead to increase of the error in some neighbouring elements. 
Because of this, the number of adaptions required by hp method per 
time step is high. Additionally, the derefinement and the p-h 
conversion steps also increase the number of grid-adaptions. It 
is clear form these figures that the p refinement method by 
itself is not able to represent the steep profile of the local 
error distribution due to large element sizes. However, when p 
-version is coupled with h - version ( by p-h conversion) the 
grid structure appears to have adequately small mesh sizes in 
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large error zones, it can also be seen from the table (4.3) that 
near the hump region, the p version results are not satisfactory 
and the local error indicator values have been polluted. But, 
after the conversion to h version and later applying the p 
version, the results in the hump region are substantially 
improved . 

An important feature of mesh distributions in the h-p 
version (Figs. 4.29 a— c ) is that for most of the time steps, the 
flame region requires not only linear degree elements but also 
higher degree elements. Also, near the flame boundaries ( where 
sharp transitions occur in the solution profile ) higher order 
elements are required. The local error indicator profiles in Figs. 
4.28 a-f also reveal this feature in the form of wide peaks 

adjacent to the sharp peaks of the flame front region ( see graphs 
for t = 0.00003, t = 0.000145, t = 0.001015 and t = 0.00431). 
The wide peaks arise due to large higher order elements; the 
narrow sharp peaks correspond to the highly refined linear 
elements . 

The effectiveness of the p-h conversion is seen in the error 
indicator for the t = 0.00011 (Fig. 4.28c ). The decrease of 

error with higher order interpolation is slow, when the number of 
elements remains same( p method). But when the p - version 
solution and mesh are converted into the equivalent h - version 
system, the fall in the error magnitude is rapid. The increase in 
the number of elements after conversion is insignificant. 
However, the refinement introduced by p-h conversion is 

non-uniform with closely spaced elements in the high gradient 
regions ( see, Figs. 4.29a-b) . Due to this, the h-p version 
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produces a well - graded mesh. The number of elements as well as 
the degree of freedom used in the h— p method are significantly 
lower than those required by the h-version, for maintaining the 
same global error criterion. It is worthy to note that the local 
error indicator values for the h— refinement procedure are normally 
smaller in magnitude than those of the p -refinement procedure. 
The reasons for this are the faster rate of convergence and the 
usage of larger number of elements in the h -method. 

In the equilibrium phase of flame movement, the h-p method 
also converges fast (within 2 adaptions ) for most of the time 
steps (see Figs. 4.28e-f and 4.29b-c). Features such as growth in 
the error band width during the ignition period (due to the growth 
in the front width ) , movement of the dense error band with 
respect to time (in order to cope up with flame movement ) etc., 
are observed similar to the h-version graphs. The mesh 
distributions also have characteristics similar to those in the 
h-method. The flame region has dense mesh of linear elements 
while higher order elements occur in the region where t.he solution 
profiles have curved shape. In some situations, (say for 
t = 0.00011), when some elements in the non reactive zone become 
coarse due to derefinement, the degree of approximation increases 
near the flame front. In the successive time steps, these coarse 
higher order elements are refined properly into a non-uniform 
distribution of linear elements ( by p to h conversion ) . In the 
flame front region, the introduction of new grid points after 
adaption have been restricted to only one in any element 
irrespective of the degree of approximation considered. This 
choice has been made in order to avoid over refinement in the 
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flame front-zone. In general , it has been observed that three steps 
of the h-p method, i. e. derefinement, increase of the 
interpolation order and p-h conversion, perform very well. In 
contrast to the h— version procedure , usually there is no movement 
of grid with adaption steps, in the h— p method. Only the degree 
of interpolation increases with adaption. Yet, the grid 
modification with time is excellent due to p— h conversion as the 
flame evolves through the ignition, equilibrium motion and 
extinction phases. 

An overview of grid evolution is shown in Figs. 4.3 0a-b, for 
the h and hp versions respectively. Further, the results also 
have been shown in the tabular form ( table 4.5 and 4.6a-b). The 
grid movement with respect to time is clearly seen to be uniform 
during the equilibrium phase. In a few places, the shift in the 
grid is not uniform due to unequal time intervals. Features such 
as the growth in flame front width during the ignition phase, 
constant front width during equilibrium movement etc. are 
observed, as discussed earlier. It is also evident that outside 
the flame region, the grid point density is slightly higher in the 
burnt, region as compared to the unburnt region. This can be 
attributed to the fact that derefinement is relatively slower than 
refinement . At every location in the burnt region, the dense 
mesh has occurred once when the flame front passed through that 
location. It takes considerable time for derefining the dense 
mesh effectively. On the other hand, in the unburnt region, the 
mesh is always coarse due to the coarse initial mesh. Comparing 
the h and hp results , it is observed that the width of the dense 
mesh region is smaller for the hp region as compared to the h 
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—version. This is explained from the fact that the edges of the 
front also require a fine mesh in h version, in order to represent 
the sharply changing temperature profile. On the other hand, for 
hp version, a few higher order elements can take care of sharp 
transition zones at the boundaries of the front. 

(III). Results of the grid speed strategy : 

The solution profiles ( temperature ) obtained by the grid 
speed strategy are compared with the h version results in Fig 
4.31a. The problem starts with a non-uniform h version mesh of 
374 elements at time t = 0.003. As discussed in the test problem 
results, the success of the grid speed strategy depends on the 
starting mesh at the initial time level. Therefore, the problem 
has to be solved by h or hp procedures until the flame front 
features are fully established and the corresponding solution and 
the mesh can be utilized for the application of the grid speed 
strategy for later times. In the figures., the shape and the 
movement of the fully established profiles for the h version and 
grid speed methods are compared over a range of time. It is 
observed that the agreement is quite satisfactory. After a 
sufficiently long application of the grid speed strategy ( with 
an optimal h -version mesh in the beginning ) collapse and 
stretching of elements occur near the boundary. This affects the 
behaviour of temperature as well as the density and therefore, 
special treatment of grid movement is done near the domain 
boundaries. As explained for the test problem, have also the nodes 
which cross the boundary are reintroduced at the beginning of the 
mesh, as can be seen from the Fig 4.31b. It is clear from this 
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figure that the nodes are distributed according to the front. 

Thus, it can be concluded for the flame problem that the grid 
speed method provides an inexpensive and efficient way of 
constructing the grids. However, its applicability is limited to 
fronts of fixed speed and shape. The flame speed predicted by the 
present grid speed strategy agrees excellently with the results of 
Bieterman and Babuska (1986) who studied the same flame 
propagation problem by h version refinement strategy. A 
comparison between our results and those of Bieterman and Babuska 
(1982 a-b) indicates that the present results are satisfactory. 
The present results have also been shown in table (4.7). 

The variation of the a posteriori error with time in the grid 
speed method is shown in the Fig 4.3 2b. Due to a good starting 
mesh and an accurate prediction of the flame speed, the a 
posteriori error remains more or less same with respect to time. 
In fact, it exhibits a small decrease in value, because of a 
slightly more optimal grid structure for the later times. 
Further, the maximum and mean errors with respect to time 
fluctuate only in the fourth decimal place and the mean of the 
errors with respect to time is almost constant. This study shows 
that the mesh constructed by the grid speed strategy is able to 
capture the front very efficiently. 

4.8.3. Computational time : 

The computational times for the h , h-p and grid speed 
algorithms have been presented in tabular form ( tables 4.8, 4.9 
and 4.10). The efficiency of these adaptive methods have been 
investigated and assessed in terms of the computer time required 






Time 

N 

Grid speed 

Apos . err 

Com. grd. 

0.003060 

374 

-141. 68312 

0.0566779 

8 .761180 

0.003160 

374 

-142 . 07661 

0.0564723 

8 .758090 

0.003260 

374 

-142 . 24686 

0.0564417 

8.754290 

0.003360 

374 

-142.34952 

0.0564358 

8.751800 

0.003460 

374 

-142 . 40625 

0.0564302 

8 .750380 

0.003560 

374 

-142 .43780 

0.0564192 

8.749590 

0.003660 

374 

-142 . 45676 

0.0564008 

8 . 749130 

0.003760 

374 

-142.46870 

0.0563758 

8.748850 

0.003860 

374 

-142.47665 

0.0563462 

8 .748690 

0.003960 

374 

-142.48224 

0.0563140 

8.748590 

0.004055 

374 

-142.48623 

0.0562824 

8.748530 

0.004160 

374 

-142 . 48964 

0.0562478 

8 .748490 

0.004260 

374 

-142.49226 

0.0562159 

8.748470 

0.004360 

374 

-142 . 49462 

0.0561335 

8.748460 

0.004465 

374 

-142.49656 

0.0561042 

8 . 748460 

0.004560 

374 

-142.49813 

0.0560837 

8.748470 

0.004660 

374 

-142.49959 

0.0560652 

8 . 748480 

0.004760 

374 

-142.50091 

0.0560485 

8 . 748490 

0.004865 

374 

-142.50216 

0.0560322 

8.748510 

0.004965 

374 

-142.50326 

0.0560172 

I 8.748520 

0.005060 

374 

-142.50424 

0.0560032 

8.748540 

0.005170 

374 

-142.50528 

0.0559871 

8.748560 

0.005270 

374 

-142.50616 

0.0559723 

8.748590 

0.005370 

374 

-142 . 50699 

0.0559574 

8 . 748610 

0.005465 

374 

-142 . 50771 

0 . 0559432 

8.748630 

0.005555 

374 

-142.50837 

0.0559295 

8.748650 


Table. 4.7. Results of the grid speed strategy 
for the flame propagation problem 


i 





223 


Dimension 

N 

Nodes 

Adp 

Non-lin 

CPU time/ 

-less time 




iter 

time step 






(in s e c o nds) 

0.000010 ' 

12 

13 


12 

0.359625 

0.000020 

23 

24 


9 

0.371920 

0.000040 

31 

32 


3 

0.230394 

0.000060 

50 

51 


8 

0.688824 

0.000100 

59 

60 


12 

1.546934 

0.000150 

141 

142 

1 

7 

2.309051 

0.000175 

232 

233 


16 

7.556550 

0.000200 

314 

315 


20 

12.734232 

0.000300 

352 

353 


18 

14.512045 

0.001100 

360 

361 


8 

6.453213 

0.001315 

356 

357 


16 

12.707541 

0.002100 

349 

350 


16 

12.491621 

0.002400 

370 

371 

1 

8 

6.581406 

0.003105 

366 

367 

2 

16 

13.257330 

0.003400 

322 

323 

1 

8 

5.780185 

0.004100 

347 

348 

1 

8 

6.229663 

0.004305 

392 

393 

2 

16 

13.908195 


Table. 4.8. CPU time of the h version algorithm 
for the flame propagation 


Dimension 

N 

Nodes 

Non-lin 

CPU time/ 

-less time 



iter 

time step 





(in seconds) 

0.003160 

374 

375 

8 

3.983766 

0.003260 

374 

375 

8 

3.985091 

0.003660 

374 

375 

8 

3 . 980361 

0.004160 

374 

375 

8 

3.986336 

0.004260 

374 

375 

8 

3.975732 

0.004660 

374 

375 

8 

3.976706 

0.005170 

374 

375 

8 

3.972914 

0.005270 

374 

375 

8 

3.988380 

: 


Table. 4.9. CPU time of the grid speed strategy 
for the flame propagation 
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Dimension 
-less time 

N 

Dof 

Adp 

Non-lin 

iter 

CPU time/ 
time step 

(in seconds) 

0.000020 

10 

16 

2 

5 

0.157503 

0.000040 

15 

26 

2 

9 

0 . 391854 

0.000060 

29 

37 

2 

6 

0.452784 

0.000100 

43 

44 

1 

6 

0.579945 

0.000150 

139 

140 

1 

7 

2.967550 

0.000175 

184 

236 

2 

11 

8 . 570375 

0.000200 

235 

236 

1 

9 

10.794509 

0.000300 

222 

282 

2 

13 

13 . 575660 

0.001200 

227 

228 

1 

9 

8 . 142457 

0.001295 

219 

220 

1 

9 

7.697207 

0.002000 

238 

239 

1 

; 

9 

8.772962 

0.002100 

234 

235 

1 

9 

8.778663 

0.003060 

218 

219 

1 

9 

7 . 018637 

0.004005 

191 

235 

2 

11 

9 . 046517 

0.004300 

200 

201 

1 

8 

6.442287 

0.005095 

199 

252 

2 

11 

10.145530 


Table. 4.10. CPU time of the h-p version algorithm 
for the flame propagation 
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to achieve a specified accuracy of the solution. 

A computational code has been written in Fortran language and 
implemented on the CONVEX C-22 0 Machine for each of the adaption 
strategies. The CPU times presented here are obtained from a 
system call routine of the Vector library available on the 
CONVEX Machine. It is a well known fact that the CPU time and the 
storage requirements are parameters which are sensitive to the 
particular implementation of a numerical algorithm. In the present 
work, a high level of code optimization has not been performed and 
therefore, the conclusions which are drawn here may be considered 
to be only qualitative. 

In the h-version algorithm, the CPU time has been calculated 
for both adaptive and non-adaptive time steps. Each of the 
important steps in the algorithm such as evaluation of elementary 
matrices, the calculation of elementary vectors during the inner 
iterations of the non-linear system of equations, decisions for 
the application of grid refinement, calculation of local and 
global error measures and the interpolation of old grid values to 
the new grid values, significantly contribute to the computational 
time. Thus, it is essential to look at the number of elements, 
number of adaptions within the particular time step and the total 
number of inner iterations in the non-linear solver, before 
interpreting the CPU time results. In table 4.8, it is seen that 
the CPU per time step is very small for initial times due to the 
small number of grid points. With grid refinement CPU increases in 
general, but occasional decrease in CPU also occurs if the number 
of adaptions and non-linear iterations to meet the convergence 
criteria are less for a particular time step( see, results for 
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t = 0.00004). A very sharp increase in CPU occurs near the 
completion of the ignition phase due to the high reaction rate 
around the hump region. Such an increase in CPU is contributed 
both by a dense refinement and a large number of adaptions and 
non— linear iterations. During the equilibrium phase of the flame 
propagation, the CPU times are more or less equal for the time 
steps without grid adaption ( e.g. for t = 0.0011, 0.0024, 0.0034 

and t = 0.0041) and similarly for the time steps with grid 

adaption ( e.g. for t = 0.001315, 0.0021, 0.003105 and 
t = 0.004305). The minor differences in the calculated CPU have a 
good correlation with the number of grid points for each time 
level . 

The CPU time results for the grid speed strategy are shown 
in Table 4.9. It is evident that the this strategy is inexpensive 
compared to the other grid refinement algorithms. The CPU time is 
nearly constant for all time steps due to fixed number of grid 
points and non-linear iterations. The minor variations in CPU can 
be attributed to the special treatments applied for the boundary 
points. 

For the hp version (see table 4. 10), the trends are quite 
similar to those of h version. However, there are two important 
points to be noted here : (i) The hp version calculations need 

more CPU for the first adaption due to the additional evaluation 
of the hierarchical basis functions and the corresponding matrix 
entries, p -h conversion steps etc. (ii) For the later adaptions, 
the CPU time increase is much less due to the small number of 
additional calculations for the hierarchical degrees of freedom. 
Another attractive feature of the h-p refinement is that the 
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number of elements required to meet the specified tolerance is 
less. 

It may be pointed here that numerical integration of 
elementary matrices is not necessary in the present work. 
However, in order to maintain the uniformity in all the adaptive 
methods, numerical integration has been adopted. The non-linear 
source term has been evaluated with eight point Gaussian 
quadrature formula. The assembled global matrix for linear degree 
basis functions is not stored for the h refinement calculations. 
Thus, optimization of the existing codes from the point of view of 
reducing storage memory and numerical computations is definitely 
possible and the CPU comparisons between these adaptive methods 
can then be made more rigorously with such optimized codes. 
However, the general features after optimization are also expected 
to be same as discussed in this section. 

4.9. Conclusions : 

In the • present chapter, the h- and hp -methods have been 
described. These error based refinement algorithms have been 
applied to the test problem and the flame problem for a stationary 
mixture The timewise variations of different error measures have 
also been analysed in detail. Apart from the h- and hp-methods, a 
dynamic adaptive strategy based on the grid velocity scheme has 
been applied. The computational results have been compared with 
analytical results for the test problem. For the flame problem a 
detailed comparative study of the h, h-p, and grid speed methods 
has been presented. 
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CHAPTER 5 

A SELF ADAPTIVE ERROR INDEX BASED REFINEMENT STRATEGY 
APPLIED TO FLAME PROPAGATION IN PRESENCE OF FLOW 

5.1. Introduction : 

This chapter describes a generalized h-ref inement algorithm 
based on a self correcting index based a posteriori error 

estimate. The adaptive algorithm has been applied to the flame 

propagation and test problems corresponding to a flowing mixture. 
Numerical results of computable error measures have been 

presented for different values of the flow velocity. The solutions 
predicted by the h- version based a-strategy are compared with the 
solutions obtained by the standard h algorithm described in 

chapter 4. 

5.2. Previous work on convective problems : Adaptive strategies : 
Convective and diffusive effects occur frequently in 

association with heat and mass transfer processes in flow 
problems. Broadly speaking, the mathematical models of these 
processes contain both a " wave like " part due to convection and 
a smooth diffusive contribution. The manner in which these terms 
are treated numerically will influence the behaviour of the 
approximate solution. It is well known that the main difficulties 
arising in the numerical solution of the convective-diffusion 
equations are due to their non-self adjoint property and the mixed 
features of elliptic, parabolic and hyperbolic character of the 
governing equations. Galerkin method, which is well suited for 
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self-adjoint problems, leads to non-physical spatial oscillations 
when applied to highly convective situations, unless relatively 
fine meshes are employed. These oscillations are more pronounced 
in the vicinity of the front where the solution gradients are 
large. Similar difficulties have been noted during the usage of 
standard finite difference procedures also, which have paved the 
way for the development of upwinding difference schemes. Thus, 
the numerical method deserves special attention while developing 
reliable adaptive schemes to convective problems. 

In order to eliminate the spurious spatial oscillations, 
various procedures have been employed over the years in finite 
difference methods, the most popular being the 'upwind ' 
differencing of the convective term. A number of numerical 
studies by finite difference methods have been reported in the 
recent past on such problems. The details of these studies have 
been summarized by Raithby, 1976; Leonard, 1981; Patel et al., 
1985). Also, several authors (Thompson, 1984; Acharya and 
Patankar, 1985, 1990; Ramos, 1985; Chargy et al . , 1990; Eriksson 

and Johnson, 1990; Rai and Anderson, 1982, Dwyer et al., 1980; 
Demkowicz and Oden, 1986 ) have proposed various adaptive 

schemes to solve convective diffusion problems. Most of these 
adaptive methods attempt to equi-distribute some measure of the 
error, but differ in the individual approaches. 

Among the various techniques employed to deal with convective 
diffusion problem, the construction of upwind finite element 
approximation via the Galerkin method with added artificial 
diffusion is an important one ( Heinrich et al., 1977; Hughes, 
Christie and Mitchell, 1978; Barrett et al., 


1978; 


1980; 
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Szymczak and Babuska , 1984; Babuska and Rheinboldt, 1980, 1982; 
Zienkiewicz and Heinrich, 1978; Carey and Oden, 1984b ). When 
upwinding is applied in finite element methods, modified weighted 
basis functions are employed to achieve the desired effect. Here, 
the test functions belong to a space different from that of 
standard shape functions. In the upwinding strategy, the 
criterion for the selection of optimal upwinding factor w has been 
either the elimination of oscillations or the attainment of exact 
nodal solutions for certain class of model problems. The concepts 
that the degree of upwinding can be continuously controlled to 
increase the accuracy of the numerical solution and that the 
oscillations can be removed by several mesh refinements, has been 
widely used in the references cited above. 

The theory developed by Babuska and Rheinboldt (1978a,b, 
1982) has also been widely used in the development of adaptive 
schemes for convective diffusive problems. These adaptive 
procedures are based on the a posteriori error estimates for a 
norm which arises from symmetrizing the bilinear form. The method 
of symmetrization has been studied by Barret and Mortan (1980) and 
has been used in the optimal FE solutions to convective diffusive 
problems. Later, Szymczak and Babuska (1984) have derived the a 
posteriori error estimates which are asymptotically exact for 
convective diffusive problems with upwinding finite elements. 
Further, the authors also suggested ways of optimising the upwind 
parameter w for certain class of problems. This a posteriori 
error estimate depends on local mesh size h and the local 
elemental error indicators calculated from the residual of the 
governing equations. The effectiveness of such an error estimate 
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has been studied by measuring' the ef f ectivity index. A variable 
local upwind ing strategy has been employed in conjunction with the 
adaptive mesh refinement by the authors Carey and Oden (1984b). 

Previous adaptive techniques ( J ohnson , 1989 ; Lohner et al . , 
1986) for convective-diffusive problems, have been based either on 
a priori interpolation error estimates for refining the mesh 
locally according to the size of the gradient or on a posteriori 
error estimates which depend only on the residue. For both the 
approaches the reliability is uncertain. In the first approach, 
the derivatives of the exact solution are used without considering 
the computed solution and in the second case, the residue may be 
large in regions of non-smoothness. Eriksson and Johnson, (1990) 
have developed adaptive streamline diffusion FE methods for time 
dependent convective diffusive problems. These adaptive methods 
employ a posteriori error estimates and the concept of space-time 
elements. The developments of Streamline Upwind Petrov-Galerkin 
method (SUPG) and the Characteristic Galerkin method (CG) have led 
to different approaches for solving convection dominated problems 
( Brooks and Hughes, 1982; and Zienkiewicz et al., 1984; Douglas 
and Russell, 1982; Oden and Demkowicz, 1986b). 

5.3. Description of the present work : 

In the past, a posteriori spatial error estimators have been 
developed for self-adjoint and non-self adjoint linear operators 
based on the error analysis for the finite dimensional approximate 
solution. However, such error analysis is generally available for 
linear operators. On a tentative basis, the same error estimates 
for linear operators have been used for non-linear problems also 
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in the past. Such usage, however is clearly ad hoc. 

In the present work, a generalized error estimator i s 
proposed which can be applied to linear or non-linear problems, it 
is only assumed that a posteriori error depends on the step size 
raised to unknown power a and the value of a is determined in the 
course of the adaption process itself. In this sense, the method 
has the feature of providing a self correcting error estimator. 
The strategy is based on computational study only and no 
mathematical theory is developed. The refinement strategy has 
been applied to the test problem and the flame propagation problem 
in the presence of mixture flow. 

In this section, we describe the test problem and the flame 
problem with convection. The details of the upwind finite element 
solution for these problems are also outlined. 

5.3.1. Description of Test problem : 

The test problem formulated here is an extension of the test 
problem described in Chapter 4 ( see eqs. 4.1a~c).. The first 
order convective term is introduced and the modified governing 
equation for the test problem can be written as 

|^u(x,t) ~( a (x)|j u(x,t)J + V ~ u ( x , t ) + b (x) u (x , t) = f(x,t) 

( x , t ) I x (0,t f ) (5.1) 
The boundary and initial conditions are same as described earlier. 
The source term f(x,t) has been chosen appropriately so that the 
exact solution of the equation (5.1.) has similar characteristics 
in the wave front region for all velocities V. For computational 
purpose, we have selected 
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a(x) = 1, b(x) =0 V X e (0,1) 

The exact solution of the eq. (5.1) is given by 

U(x,t) = 1/2 + 1/2 tanh [2/3 (x-Ct) ] 

with same values of /3 and C as described in the Chapter 4. (see 

eq. (4.2) ) . 


5.3.2. Description of Flame problem in presence of flow : 

The governing equations of the flame propagation problem in 
the presence of flow can be written as follows : 


d T 

a t 

a_p 
a t 


a 2 T 
a x 2 

a 2 p 


v 


V 


a t 
a x 

a p 
a x 


= f(P,T) 


= - f(P,T) J 


t > 0, x e (0,1) 


ax ax (5.2) 

where V is the prescribed fluid velocity. The boundary and initial 
conditions for this problem have already been defined in equations 
(3.1b-f) of Chapter 3. The above equations can be expressed in 
the more general form as follows : 

(u i J t “ k(x)( Ui ) ] + V (u ) = fj (x,t,u) ; (x,t) efix (0,t f ] 

> X 1 = 1, MPDE 


(5.3) 


During the computations, a t (x) has been set equal to unity in all 
the equations, for the sake of simplicity. 


5.3.3. Basic formulation: 

This section describes the mathematical frame work in which 
the above test problem (eq. 5.1) and flame problem (eg* 5.2) will 
be studied. The notation used and the definitions of the 
variables employed here, are analogous to those defined in 
Chapters 3 and 4. 
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Semi discrete approximation : The finite element formulation and 
the method of solution have been explained in detail Chapter 3. 
Here, only the important steps in the FE formulation have been 
discussed. Following precisely the same procedures employed in 
Chapter 3, the semi discrete approximation for the eqs. (5.3) can 
be written as : 


< : 


5_U (.,t) 

at 


Vj> +< a 


1 


£U 

ax 


a v 
ax 


1 > +< v 


au (x , t) v 

3x ' 1 


> 


= < fj (*,t,U) ,v x > 

for all U -V, € s[ ,k •, t e (0,t ], and 1 = i...,npde. 

1 1 h f 

(5.4a) 

with suitably chosen initial conditions. Expressing U (x,t) in 
terms of basis, the above system of equations can also be 
written in a matrix form as 



where 



" 


| 

CV’ 

1 I 


(uct)} = |f i (U [ t ] ) 


- 

for 1 ■ 1 , NPDE 


(5.4b) 


[ K l] = an<j [ K f| = * \). (1 * j * v 

are the matrices with entries 

(kJ) = <aV\p , ViJj > 

1 ' 1 j 1 *i 1 ' r 1 j 

( K f) . , = <a , \fj >. 

1 1J 1 ^li' ij 


and 


k 

Here, is the dimension of the finite dimensional space S^’ 

and UJt] ( 1 < npde ) is the vector of nodal degrees of 

freedom which is determined as the solution of the system of 
time dependent differential equations (5.4b). The other terms in 
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the above equation are defined in the same fashion as in 
Chapter 3 . 

upwinding formulation ; 

In the present study , the concept of upwinding finite element 
technique has been employed in the solution procedure. A linear 
combination of special test functions (upwinding basis functions ) 
with piecewise linear basis functions have been used for test 
space of solutions in the weak formulation (see Szymczak and 
Babuska, 1984) . This is achieved by adding a quadratic term, 
multiplied by some parameter w to each linear basis function of 
the test space. The construction of upwind finite basis functions 
for a particular element is explained below : 

Let % (?) and ^ (?) ; -1 ^ ^ 1; i = i, npde, be the 

1 > 1 i , c. 

usual piece-wise linear basis functions on the master element 

Q n = [-1,1]. Using standard transformation procedures, these local 

basis functions can be used to obtain the global basis functions 

^ (x) and i jj (x). The upwind weighted element functions (test 

1 » 1 1,2 

functions for the element Q ) can be written as follows. 

n 

*ij?) = + w \ 2 (?) 

\ a <5> = \ 2 W ~ w 

(5.5) 

where w is a scaling factor that specifies the amount of upwind 
bias desired. Here, the product ¥ ^ x (?) ^ 12 (?) i s the additional 
quadratic form in the test function and it vanishes at the end 
points of fi (Fig. 5.1). Further, the sign in the expression 

n 

indicates the upwind direction which may vary according to the 
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problem. For the flame problem in which the front propagates from 
right to left, the upwind direction is defined from right to the 
left. Thus, at every node, a quadratic bias is added or subtracted 
from the standard linear basis function corresponding to the 
region ahead or behind the node respectively, in the direction of 
propagation. Using the upwind biased basis functions ( i (?) 

\ 1 > 1 vs Jr 

-(C)) as test functions for the master element Q , one can 
construct the global test functions with the desired properties 
(see Fig. 5 . 1 .) . 

Representation of the solution space : 

In order to obtain a finite dimensional approximate solution 
of equation ( 5 . 3 ) , the finite dimensional solution spaces for the 
test and trial functions must be specified. For the sake of 
convenience, we recall the definition of S^’^(A) which has been 
used in the Chapter 4 . 


• r k . # 

Trial space : We denote S, ’ (A) as the trial space which consists 

— h 

of all (r- 1 ) times continually differentiable functions on Q 

for which the restriction to Q ,j = 0,1,2 ,. .jV - 1 is a polynomial 

j h 

of degree k -1 and k a 2 r, j = 0,1,2 1. 


1c * 

Test space : The test space S ’ (A) is spanned by two sets of 

basis functions, namely, the linear basis functions 



and the upwinding basis functions 



where, the upwinding basis functions 



j = 1,2 N 

n 

# (x) are 

1, j v 


constructed by 
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the procedure described above. It is to be noted that the vector 
£ in the test space and k in the trial space are not equal, as 
their components may differ from each other due to different 
degrees of approximation considered locally in each element of the 
solution domain. In general, the basis functions in the test 
space are denoted by 

*, - { : J = '•* # h 

where x (x) is a linear combination of 

i,J v ' 

V. ,(x) and $ (x) 

(5.6) 

Implementation of Upwinding : 

To generate the desired upwind model, the test space 

r 

S’ (A) described above is considered for the purpose of 
constructing the approximation. Using the following test functions, 
V = x ; x eS r ’k(A), for i = npde 

1 1 1 h 

the semi discrete equations in the upwind formulation can be 
written as : 


3 U 
3 


t 


( • ft) 

f 


+< 


au (x,t) 
i ax 


3_x 

ax 


i >+< v 


3 U ( x , t ) a 
ax ' ? 


= < f j (x, t,u) ,<c l > 

for all U € S r ’k(A) ; 1 = 1 ,... , N P D E . 

1 h 

(5.7a) 

with suitably chosen initial conditions. The resulting 
equations (5.7a) can be written in matrix form as : 
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“ "I Fj (U[t] ) j- + w |f i (U[t])J- 

for 1 = 1 , 2 ,... NPDE 

( 5 . 7 b) 

In the above equation, the quantities with have been 

contributed by the product terms of the upwind basis functions 

while the other matrix contributions are due to standard basis 
functions in the test space. Also, the superscripts i and 2 
indicate diffusive and convective contributions respectively. 

The upwind matrices are given by 

[*!] - [«?] - 

(1 =s i < N ) , (1 £ j s N ) 

n n 

with entries 

(K 1 ) = <a Vi/f . Vm > 

v rij 1 *11' 1 j 

(K 2 ) = <a Vi/f , <c >. 

' rij 1 Y u' 1 j 

and the upwind vector (U[t] ) j- has entries 

= <f i (x,t,U[t] , <c>. 

5 . 4 . Description of present a posteriori error estimate : 

A major contribution of the present work relates to the 

development of a generalized equi-distribution strategy 
(designated as the a-strategy ) using for local error indicator 
and the a posteriori error. The error is estimated in terms of the 
residual of the finite element solution and the local mesh size 
parameter raised to an unknown index a. The only assumption 
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invoked for the application of the a-strategy is that the error 
estimator depends on the product of a measure of the local residue 
in some properly chosen norm and a positive power of the local 
mesh size. Based on these requirements, the a posteriori error 
estimate g i,m (t,a) at t = t for the m th adaption can be written 
as follows : 


S i,m (t,a) = 


N 


i , m 


n = 1 


1/2 kf 1 , m 

f e NPDE 


■ II 


n 1 / 2 


i % m m | % 

V. (t,oc) 

1 , n 


n = X 1=1 


(5.8a) 


where the error indicator V ’" 1 (t) is given by 

1 , n 


l,n 



n 


( r !’n (t ' S) ^ dS 
l,n 

n 1 ’ m 


(5.8b) 


Here, r|’ m (t,x) is the residue of the finite element solution in 
n th element which is computable for each t > 0 on the sub 
interval (x , x ) . 

v n ' n+ 1 ' 


5.5. Scope of the present study : 

The aim is to develop a computational procedure which 
meets the following goals : 

(Gl) . To develop a versatile, error based refinement strategy 
which has the capability of producing sufficiently accurate 
numerical solutions of non-linear parabolic partial 
differential equations. The accuracy is to be measured in 
some properly chosen norm associated with the problem. 

(G2) . To develop FE solutions for the flame and test problems in 
the presence of flow which satisfy a prescribed a posteriori 
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error criterion. 

(G3 ) . To obtain the proper distribution of grid points for 
h— version oc strategy at each instant of time for the 
physical problems stated in (G2) . 

(G4) . To present a detailed study of behaviour of various error 
measures for this new algorithm and compare the relative 
performance of the proposed error estimate with that of 
standard h refinement strategy. 


5.5.1. Mathematical statement : 


In this 

subsection, the 

objectives 

are expressed in 

the 

mathematical 

frame 

work. 

Firstly, the 

global 

residue 

norm 

SR i,m ( t) is expressed 

in terms 

of elemental 

residue 

norms 5R 1 ’ 

n 

m (t) 


as 


* l ’ m 


N 


i , m 


NPDE 


fR i,m (t) 


n 


II 


n = 1 


i . m f i » 2 j 

r ’ (s,t) ds 

l,n 


1=1 n = 1 Q i,m 


We define the relative global residue ft*’"'(t) 

Rl L 

adaption at time t = t ; as 


(5.9a) 

for m th ( m > l) 




i , m 
REL 


(t) = 


ft i,m (t) - K l,B-1 (t)j / K i,m (t) 


(5.9b) 


The immediate goal is achieved in two steps. Given a problem (P) 
with computed solution U(x,t) and certain measure of error the 
tasks to be completed are : 


(a). For the initial time interval (tg,^], design an adaptive 
algorithm based on equi-distribution of local error 
indicator for constructing a finite element mesh ' A' such 
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that the residues ft 


1 , m-1 


, i , m 


(t) and ft ' (t) for two successive 


adaptions (m-i , m) satisfy the relation 


ft 1 ’ m (t) 5 

R E L ' * 


REL 


and determine the index a by setting 


S 


i , m 


(t,Q£) “ 5,= 0 


(5.ioa) 


(5.10b) 


where g 1,m (t,cc) is defined in the equation (5.8b). Further, 
the equi-distribution is performed with the help of the 
statistical procedure. 

(b). For the subsequent time intervals, the task is to design an 
adaptive algorithm for constructing the finite element mesh 
'A' such that a posteriori error S i,m (t,a) satisfies the 
relation 

i , m 


Sj- ’“'(t,a) £ 5 2 V t e (t t ,t f ) 


(5.11a) 

where a is known apriorily or determined from the relation 
£ i,m (t,a) - 6 i = 0 (5.11b) 

The decisions concerning mesh modification are aimed at 
achieving (5.11a) at different time instants .Whenever 
g i,m (t,a) > 5 , the grid adaption is carried • out only once 

i in 

and at the same time & ’ (t,a) is reset to 5 i for the 
determination of a on the refined mesh. Here, S and 5 are 
the prescribed lower and upper tolerance parameters for 
bounding the a posteriori error. 


The above aim has the following detailed description for each 
time interval (t jf t i+i ] of total time interval [t ,t f ). 

The adaptive algorithm seeks to construct a mesh A I,M with mesh 
size vector h l ’ M and the corresponding FE discrete solution 
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( x >t); (l<l< NPDE) 


such that 


i , M 


5 - g (t,a) s 5 or x ’ M (t) s 5 

2 REL' ' REL 

for t € (t i# t i+i ] (5.12a) 

where a is known apriorily or determined from the relation for 
all adaptions m = l, m 

gi ’ m (t.a) - 5 j = 0 (5.12b) 

This is achieved by constructing a sequence of meshes 

A 1 ’ 1 A 1 * 2 A 1 ’ 3 ,i,« 

A ' A F A / /A (5.12c) 

of mesh size vector 


h 1 ’ 1 , h i,z , h i>3 , 


/ ™ 


with corresponding computed solutions 


(5 . 12d) 


uj’ 1 (x,t) , uj’ 2 (x,t) , 


, Uj lM (X,t); 


u 1 ’ m (x,t) € y i,m ( A) 


where A 1 ’ m + 1 is constructed from 

U i,m (x,t) e y i,m (A) 

by equi-distribution of the elemental error 
contributions 

E i,m (t,a) ; (i s n s A' i,m ) 

n e 

to the a posteriori error 

£ i,m (t,a) , t e (t jf t. +i J . 


(5. 12e) 


(5 . 12f ) 
indicator 


Here, the solution U i,m (x,t) is obtained by solving 


the 



244 


semidiscrete problem 


a u; ,m (x,t) i, 

< 8t 1 ' l 


> + <a 


l d x 


Uj’ m (x,t) 


a x 

sTx 


i , m 

i > 


+ v < 


a x 


U, ,m (x,t) , > = < f i (x / t,U i,m (x,t) ) ,a:|’ m (x,t)> 


V u|' m (x,t) e y i,m (t) V a; ] i,m (x,t) e if i,m (t) 

(1<i<npde ), and m = 

(x,t) e Q x (t t ,t i+i ] 


( 5. I2g) 


with appropriate boundary and initial conditions. In the equation 

(5.12a), the stopping criterion is based on residue for initial 

times t 3 t i ] and on a posteriori error for later times 

(t < t s t ] . 
v i f J 


5.6. h-version a strategy : 

In this section, the a strategy based on h-version grid 
refinement is presented. Various steps involved in the finite 
element solution and the mesh modification are described in 
detail . 


5.6.1. Description of the algorithm : 

The timewise evolution of the FE solution and the mesh 
distribution with time or adaption level are described here. At 
the beginning of each time step, the initial data has been 
specified as described in the h version algorithm of Chapter 4. 
After each adaption , the mesh distribution and the starting 
guess solution change, while the initial condition requires 
interpolation to the new mesh. On the basis of the computed 
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Fig. 5.2. Flow chart for h-version based a strategy 
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finite element solution, a posteriori error ( or the relative 
global residue for the initial time step ) is estimated with the 
help of local elemental error indicators . The decisions for mesh 
modification are then generated by the equi-distribution of the 
local error indicator E (t,oc) . For the initial time interval the 

n 

adaptive grid is applied till two successive adaptions satisfy the 
tolerance criterion of eq. 5.10a. Here for each adaption the value 
of a is calculated and the smallest value of a is used in the 
error indicator expression for the mesh refinement decisions. 

At later times, whenever the a posteriori error exceeds the 
upper tolerance parameter ( 8^) grid adaption begins, using the 
currently available value of the power index a for evaluating the 
error indicators through equation (5.8b). The grid adaption is 
done only once. On the refined mesh the finite element solution is 
calculated and the new value of a is evaluated from the eq. 
(5.11b) . Further, at this stage, the value of a .is updated by 
choosing the the smallest value from the previously available and 
the newly calculated values. This updated a is used while 
calculating the a posteriori error g(t,a) at t = t j+i - A l so / the 
posteriori error is reset equal to to lower tolerance parameter 
(6 i ) . The solution at this stage is accepted as the computed 
solution for the updated time level. At later times, grid adaption 
is not carried out whenever the a posteriori error lies between 
two tolerance limits. The sequence of steps performed for each 
adaption cycle are enumerated below. These steps are also shown in 
the Fig. 5.2, in the form of flow chart. 
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Step 1 : Calculate the finite element solution at t = t 

i 

Calculate the FE solution at t - t. tor m th adaption on 
the mesh A 1 ’" with given Initial condition and 
starting guess solution. 

Step 2 : Check for nature of adaption for given time interval 

If t e ( t Q , t^ then go to step 3a for the evaluation of 
relative global residue *£“(t) and the index a. 

If t £ ( t Q , t i J then go to step 3b for the evaluation of 
a posteriori error & ’ (t,a) and determination of a. 

Step 3a : Calculate R^“(t), index a and E*’ m (t,a), (i s n s ^ i>m ) 
Calculate •^■ R ’ ZL (t) from equation (5.9b), compute 

index a from equation (5.10b) and elemental error 
indicators E^’ m (t,a) and go to step 4a. (For * = i, 
choose (t) = 0 ) 

Step 3b : Calculate index a and E i,m (t,o:), (i < n < ^ i,rn ) 

n e 

If m = 1 (starting adaption) then 

Calculate the elemental errors E 1,m (t,a), (is n=;A' l ’ m ) 

n e 

Calculate the a posteriori error estimate <s 1,m (t,a) 
with updated a and go to step 4b. 

If m = 2 (final adaption) then 

solve for new a from the eq. g i,m (t,a) -5 =0, 

reset the a posteriori error g i,ni (t,a) to 5 . 

Select the least of new and previous values of a 
and go to step 7 for time marching. 

Step 4 : Check for the application of adaptive grid . 

(a) : if t e (t ,t ] and ft i ’ ra (t) s s then stop the 
adaption and go to step 7 for time marching, 
if t e (t ,t ] and ft*’ m (t) > 8 then go to step 5 

01 REL REL 
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for adaption. 

(b) : if t s (t i> t f _] and 5 ^ g 1 ’ m (t,a) ^ 5 2 then 
go to step 7 for time marching. 
if t e (t ,t ] and & 1 ’ m (t,a) > S 2 then 
go to step 5 for adaption. 

Step 5. : Implementation of adaptive h-refinement scheme 


For 

for any 

y-v i. y m a i « in * /-» . *pi i > hi / * % 

fi ’ € A , if MINTOL < E ( t , a ) 

n n 

< / l i ’ m 

then 

RETAIN 

the element as it is. Here, n 

IS 

the 

mean value of the error 

indicators E 1,m 

n 

(t,a) 

and 

MINTOL is 

the specified 

tolerance for 

node 

removal . mintol 

is dynamically 

set as 



MINTOL 

d *n l,m where (o < 

d < 1) . 



(ii) . For any pair of elements , (Q*’ m ,Q^™) if 

E i ’ m (t,a) A E i,m (t,a)^ mintol then REMOVE the 

n n+ 1 

common node between these elements. 

(iii) .For for any Q i,m e A 1 ’™, if E 1,rn (t,a) > n. i,m then 

n n 

REFINE the element according to statistical 

procedure. 

(iv) . Perform the above steps 5 (i) , (ii) , (iii) and (iv) 

for all elements on the mesh . 

(v) . Generate initial data on the new mesh A 1,m+1 

increase the adaption counter (m) by one, and go to 
step 1. 

Step 6 : Time marching : 

if t < t f then set t j+i = t } + At; construct and 
generate the initial data on the coarse mesh A 1 + 1 ’ 1 at 
t = t l+1 > reset counter (m) to 1 and go to step 1 
if t = t f then stop the algorithm. 
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5.6.2. Important features of Algorithm : 

The procedure for evaluating the initial condition, starting 
guess solution at any adaption and time level are similar to 
those used in the h -refinement algorithm of Chapter 4. Hence the 
details of these procedures are not discussed here. 

Finite element solution : 

The matrix contributions [M] , [K 1 ] and [K 2 ] in the equation 

(5.7b) are due to the additional terms in the test functions with 
quadratic bias. These modifications arising from upwinding 
formulation are easily implemented at the element level and the 
resulting matrix system (5.7b) is solved by the iterative method 
described in Chapter 3 . 

Decisions for mesh refinement : 

In the following section, the adaptive mesh modification 
based on the newly proposed a posteriori error estimate and 
some other important aspects of the adaptive mesh generation 
procedure are discussed. 

For the initial time interval, the adaptive grid refinement 
is carried out based on the equi-distribution of local error 
indicators, till the global residue for two successive adaptions 
satisfies the eq. 5.10a. After a sufficient number of adaptions 
the global residue and the index oc attain stable values. 

For s ub sequent time intervals, when the a posteriori error 
exceeds the upper tolerance & 2 , mesh refinement is performed 
according to statistical procedure and the value of ec is updated 
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on the refined mesh using equation (5.11b) . Also, the a posteriori 
error is set equal to the lower tolerance limit <5 ^ . Grid adaption 
is not carried out whenever the a posteriori error lies between 
two bounds. In actual computations, at later times, the number 
of adaption steps is always restricted to two, since the a 
posteriori error turns out to be nearly equal to the lower 

tolerance parameter after one adaption. The evolution of a 
posteriori error with adaption and time level is shown 
schematically in Fig. 5.3. 

In the de-refinement process, the criterion for the removal 
of nodes is set dynamically as a given fraction of the mean of 
local error indicators. A constant 'd' ( 0 < d < 1) has been 

chosen such that the criterion for the removal of nodes is 

given by mintol = d*n. i,ra . A low value of the constant will 

result in slow removal of nodes in comparison to the stuffing of 
nodes. In the computations of the present study the value of d has 
been chosen as 0.75. 


Calculation of index a. : 

Let us suppose that the adaptive grid has been applied on the 
mesh A i,m at t = t^or the m th adaption. We calculate the index 
a from the equation 


S l ’ m (t,a) -5 = 0 


N 


i , m 

e NPDE 


i*e. , 


z i h “ 


n = 1 1=1 


1 % m , . . 2 •« 

r (s,t) ds 

l,n 


5 = 0 = V(a) (say) 


Q 


1 , m 


(5.13) 

The root oc is calculated by the bisection method. Since the step 
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Results have been predicted for various values of flow velocity. 
The situation without flow ( V = 0) has also been considered as a 
special case, for the sake of comparison with the results 
presented in Chapter 4. The nature of the solution profiles, 
behaviour of the global and local error measures and the variation 
of the exponent a are discussed. 

5.7.1. Application of a strategy to the test problem : 

(a). Solution profiles and grid distributions : 

The exact and the computed solutions of the test problem are 
compared at various time levels in Fig. 5.4a-c. for three 
different flow velocities. The same exact solutions have been 
used for all these velocities, by selecting the source term f(x,t) 
of eqn. (5.1) appropriately. It is observed that the exact and 
the computed solutions match excellently at each time, during the 
initial transient stage as well as the uniform movement phase of 
the front. For the first few time steps , in order to compensate 
for the inaccuracies arising from the coarse mesh, the initial 
condition for each time step has been taken from the corresponding 
exact solutions. 

The front shape and its movement are predicted very well by 
the a-strategy analogous to the refinement strategies based on the 
theoretical error estimates of Bieterman and Babuska (1982a, b, 
1986) described in Chapter 4. However, it should be noted that 
these theoretical error estimates are not applicable to the 
non-linear non self adjoint operator, of the present problem. In 
this sense, the a-strategy could be a powerful tool for generating 
meshes in an adaptive manner for any general problem. The tables 
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associated with Fig. 5.4a-c show the number of elements at 
various time levels. At time t = 0, a uniform mesh of 10 elements 
has been used in each case. It is seen that the number of 
elements increases rapidly during the initial transient stage and 
rather slowly when the front has a uniform movement phase. These 
features of grid refinement are closely related to the improved 
estimation of the exponent a with each grid adaption. Indeed the 
a-estimate improves rapidly during the initial adaptions and 
slowly in the later stages. More discussion on these trends of 
the exponent is presented subsequently. During the initial time 
interval, the computed solution is accepted if the relative 
residue satisfies the tolerance ( 10~ 2 ) criterion. 

In general, the number of elements at any time is somewhat 
less for the a-strategy as compared to the procedures discussed in 
Chapter 4. With respect to the flow velocity variation, the 
number of elements is more or less the same for low velocity 
values. For high flow velocity values (V = 10.0), however, the 
number of elements is significantly more. The grid distributions 
for various flow velocity values and time instants are shown in 
Fig. 5.5a-c. It is evident from the figure that the cc-strategy 
provides dense refinement in the front zone and that the grid 
distribution faithfully follows the front movement with respect to 
time. It appears, however, that the magnitude of error is 

different for each flow velocity; this may be the reason why the 
number of elements is not the same, although the same exact 
solution has been used in all the cases. The plots also indicate 
that grid refinement occurs in anticipation of front arrival in 
regions ahead of the front, while gradual derefinement occurs in 
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the regions behind the front. These features are highly essential 
for any adaptive strategy which is to be used for the capturing of 
mobile high gradient fronts. 

The adaptive algorithm automatically removes, retains and 
stuffs nodes simultaneously. It is however, possible that some 
unnecessary nodes exist in regions where no solution activity is 
taking place, as observed from the Figs. 5.6a-b. Also, at few time 
instants, refinement is slightly less dense, even in the flame 
front region which can be seen from the figures. Except for some 
minor exceptions, the figures clearly illustrate the ability of 
the adaptive scheme to concentrate grid points around the 
characteristic wave, where all the activity takes place. The side 
tables provide the information of refinement and de-refinement of 
elements in the adaption. 

In the present work, removal of nodes has been carried out if 
the error of the neighbouring elements goes below a prescribed 
fraction of the mean error. Thus, the minimum tolerance limit for 
removal of nodes in the passive regions, is set dynamically. The 
removal of nodes in the de-refinement process is very fast if the 
fraction of mean error (d) is chosen appropriately. The value of 
fraction d has been taken as 0.75. Due to this relatively large 
cut off value, however, mesh modification is not done properly at 
some time instants, although the trend is alright and high 
gradient regions are refined as expected (Fig. 5.6a). On the other 
hand, over refinement of the grid occurs some times due to the 
dense stuffing of nodes by the statistical refinement procedure. 
However, this can easily be handled by specifying a minimum grid 
spacing. Further, for the first few adaptions, the resolution of 
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front region is not proper because of coarse initial grid. The 
main reason for this is the improper evaluation of the index a 
which in turn significantly influences the magnitude of the local 
error indicator. However, this problem is eliminated in the later 
adaptions as the grid becomes more refined, (see Fig. 5.6a-b). At 
later stages of grid adaption, the mean error is very small and it 
is prone to be affected by round off error. Thus, at few 
instants, elements which need not be refined are also divided 
further. In the computations for the a strategy, a large time step 
value ( At = 10 4 ) has been used. Because of this, the local 
elemental errors do not behave perfectly well, especially adjacent 
to the moving front. Large values of At, cut-off fraction for 
defining mintol etc. have been purposefully used here, to subject 
the a strategy to a more stringent test. If these values are 
lowered suitably, excellent results can be obtained as in Chapter 
4, at extra computational cost. It can be seen that the number of 
elements required to satisfy a similar error tolerance criterion 
as in Chapter 4, is much less for the a strategy. Figs. 5.6a-b, it 
is clear that the nodes are distributed in anticipation of the 
front motion. 

(b). Behaviour of index a : 

In Fig. 5.7, the exponent of a has been plotted against time, 
for three different flow velocities. The timewise evolution of a 
including the variations over each adaption step of the time 
increments At, is also depicted. It is observed that a has a high 
value in the beginning when the global residue (error) is large 
and the mesh is coarse. For such coarse initial meshes, a good 
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value of a is obtained by using an adaption stopping criterion 
based on global relative residue. This is equivalent to reducing 
the error in the solution at initial time. with successive 
adaptions, the value of a decreases rapidly and attains a stable 
value for refined meshes. This process can be likened to the 
approach towards the asymptotically valid error estimates of 

Babuska and Rheinboldt ( 1982a, b) for h » o. Indeed for refined 

meshes, cc approaches the value of 2, as in the error estimator 
expressions described in Chapter 4. This proves the self 
correcting feature of the a-strategy, which leads to improved 
estimation of the computable error after each grid adaption cycle, 
until the computable error approaches the true error for the 
problem. . The minor oscillations seen in the a value with respect 
to adaption level within each time step, arise due to the front 
movement during the time increment At. Also, the choice of a large 
time ( At =* 10~ 4 ) step may be another reason for these 
oscillations. In fact, the error tends to increase as the front 
moves out of an optimal grid distribution at a particular time 
instant. This increase in error tends to increase the a value 
slightly. However, in the next adaption cycle, the grid 
distribution is adjusted again to the new front location and the 
value of a is corrected back to more or less the original value. 
The overall variation of a with respect to time follows a 
monotonical ly decreasing trend, confirming the progressively 
improved error estimation feature of the present strategy. For 
all the flow velocities considered, the behavior of a essentially 
remains the same. These trends confirm the proposal of Szymczak 
and Babuska (1984) that in the asymptotic limit of h > 0, the 
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computable error is well estimated by their error estimator 
expression for the linear convective-diffusive problems. 

(c). Behaviour of residue and error estimates : 

In Fig. 5.8, the variations of the a posteriori error and the 
gradient error with time are shown. Since two limits (5 j and 8 2 ) 
have been employed for starting and stopping the grid adaption, 
the a posteriori error oscillates with respect to time within 
these two bounds. Occasionally, during the course of a particular 
adaption or time step, the error marginally exceeds the upper 
limit or dips below the lower limit. The plots of gradient error 
indicate that for refined meshes, maintaining the a posteriori 
error within certain bounds implies the limiting of the gradient 
error as well. Moreover, as time proceeds, the mean value of the 
error decreases and eventually stabilizes. 

In Fig. 5.9, the variations of the maximum and the mean local 
error indicators are plotted with time. The instantaneous values 
of these error indicators oscillate with respect to time, as a 
consequence of limiting the a posteriori error within two bounds. 
However, their average values decrease with time monotonically 
indicating successive improvement in the computed error estimate, 
the associated grid distribution and the predicted solution. The 
magnitudes of the maximum and the mean error indicators are almost 
of the same order, implying that the statistical procedure is very 
effective in equi-distributing the local error. With respect to 
flow velocity, the computed error values show some variation, 
although the range of minimum does not differ much. 

The maximum and the mean values of the local elemental 
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residues are plotted against time in Fig. 5 . 10 . These elemental 
residues also decrease rapidly in the beginning approaching 
towards stable values for refined meshes, similar to the features 
exhibited by the other error measures. 

The variation of the global residue and the global solution 
gradient (in the respective norms) are shown in Fig. 5 . 11 . For 
initial times , these norms increase due to the growth of the 
front. In the later stage, when the front is fully developed and 
its velocity is constant, the norms of the residue and the 
gradient also attain constant values. 

The spatial distributions of the local error indicators are 
plotted in Figs. 5.12-5.14 for a flow velocity value of V = 5.0. 
These figures show that the magnitude of the local error indicator 
decreases with grid adaption, as expected. The number of 
adaptions required to meet the a posteriori error tolerance 
criterion is large in the beginning, due to the coarse initial 
mesh. For later times, however, just two adaptions are 
sufficient. The band of sharp error peaks (representing the zone 
which contributes primarily to the global error) moves with time, 
in response to a similar movement in the front location. A few 
error peaks of large width are also seen in the figures, at some 
time instants. These imply that one or two elements lying outside 
the front region have not been refined enough; but the 
contribution of such elements to the global error is negligible 
and it does not matter even if the local error has not been 
equi-distributed well over these elements. 

(d). Behaviour of error measures for the case V - 0 : 
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The exact and the computed solutions of the test problem are 
compared at various time levels in Fig. 5.15 for the case V = 0. 
It is observed that both the solutions match excellently at each 
time, during the front movement indicating that the refinement 
process works very well and is very accurate. The grid 
distributions for various time instants are shown in Fig. 5.16. 
It is evident that the a— strategy provides dense refinement in the 
front zone and that the grid distribution faithfully follows the 
front movement with respect to time. The other features are 
similar to those observed for the non-zero velocity case. The 
spatial distribution of local error for the non-convective 
situation (V = 0) are shown in the Figs. 5.17-5.19. Here also, 
the band of sharp error peaks moves with time, in response to a 
similar movement in the front location, as observed in the 
convective case. The mean and the maximum of the residuals as well 
as the error indicators have been shown in the Fig. 5.20. These 
also have similar characteristics as discussed in the convective 
problem. Also, the figure shows that the value of index a as 
expected, is nearly two during the uniform front motion. Further, 
the results for the case V = 0.0 are presented in the table 5.1. 

5.7.2. Application of a-strategy to the flame propagation problem : 

For the flame propagation study, some results have been 
predicted using uniform meshes in addition to the adaptive grids 
of the h-version a-strategy. This is in view of the fact that 
analytical solutions are not available for the non-linear flame 
propagation problem and uniform mesh solutions (on sufficiently 
refined grids) can be used for the sake of comparison. Moreover, 
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they are useful for analyzing the effectiveness of grid adaption. 
In the ensuing sections, the behavior of the solution profiles, 
grid distributions and various error measures are presented and 
discussed. 

(a). Solution profiles and grid distributions : 

In Fig. 5.21, the solution profiles of a - strategy are 
compared with those obtained using uniform meshes. The upwinding 
finite element procedure has been used for here to obtain the 
numerical solutions. During the ignition phase, even a coarse mesh 
(of just 18 elements) used in conjunction with the a-strategy 
provides as accurate a result as the uniform mesh of 400 elements 
for V = 5.0. When the flame front moves at a constant speed, the 
predictions over the uniform mesh are not accurate enough. For a 
200 element uniform mesh, the results are very inaccurate, 
especially in estimating the correct propagation speed of the 
front. The shape of the front is predicted reasonably well. The 
results of a 400 element uniform mesh are closer to th'ose obtained 
by the a-strategy; even then small inaccuracies still remain in 
the predictions of the uniform mesh. This is in view of the fact 
that grid point locations are optimal at each time instant in the 
error based a-strategy, while in the uniform mesh many elements 
are located in regions where they are not needed. More discussion 
on the high accuracy of predictions for the a-strategy is provided 
subsequently in connection with Fig. 5.23. Thus the adaptive 
strategy provides a superior performance in terms of good solution 
accuracy, even with less number of elements. It must, however, be 
mentioned here that the number of elements increases with time for 
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the a-strategy till the error estimate asymptotically stabilizes; 
the total number of elements required to meet the a posteriori 
error criterion after the stabilization of a, depends upon the 
values of the error limits (5 and 8 2 ) . With the proper choice of 
these limits, quite accurate results can be obtained with moderate 
computational effort. 

In order to test the general applicability of a-strategy , an 
extreme case of V — 200.0 has been considered in Fig. 5.22. Due 
to the very high flow velocity, it is expected that the front 
motion will be completely arrested in this situation. Indeed for 
large time, steady profiles of steep variation in the temperature 
and density are observed, adjacent to x = 1.0, indicating attached 
flame at this boundary. 

The solution profiles obtained by the a-strategy are compared 
with those predicted by the conventional h-version method for the 
no flow situation of V = 0 in Fig. 5.23. Excellent agreement is 
seen in Fig. 5.23 and any minor deviations present rare solely 
attributable to the inaccuracies of the uniform mesh solutions. 
The number of elements is some what higher for the h-version 
method in the beginning. However, this trend reverses with 
a-strategy requiring more elements, for later times. Initially, 
on a coarse mesh, the value of a is high and it stabilizes slowly 
with time. Due to this, the a-strategy employs less elements (or 
larger average step size ) to meet the error criterion in the 
beginning. However, when the value of a is stabilized, the 
a-strategy needs more elements, due to the fact that it uses two 
error limits (5 j = 0.50, and 5 2 = 0.75), with the upper limit S 2 
coinciding with the error tolerance of the h-version method (S 2 
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= 0.075). Since it is required that a posteriori error should 

lie between 8 and 5 (with 5 < 8 ), the application of the error 
1 2 12 

tolerance is more stringent in the a-strategy, which explains the 
larger number of elements during the later stage of front 
propagation. 

In Figs. 5.24 and 5.25, the grid distributions are shown for 
the velocity values of V = 0.0 (see Table 5.2) and V = 5.0 at 
various times. All the features of the adaptive meshes such as 
dense refinement in the front zone, movement of the grid points 
along with the front etc., are clearly observed. Indeed the fact 
that the value of the exponent a is not initially stabilized does 
not seem to affect the nature of mesh refinement. Thus there is 
a strong case for using the a-strategy in nonlinear problems for 
which theoretical estimates of error are not available and the a 
posteriori improvement of error estimation may be the only 
suitable option. Comparing the convective and non-convective 
situations it can be said that since the flow velocity direction 
is opposite to the front propagation, the. effective flame speed is 
reduced when the flow is present; therefore less number of grid 
points are sufficient to meet the error criterion. This feature is 
clearly evident for an extremely high flow velocity V = 2 00.0 as 
seen from the table provided with Fig. 5.22. It appears that one 
of the main reasons why a large number of grid points are 
required to study flame propagation is the high speed of the flame 
front, due to which a zone thicker than the flame front has to be 
refined. Therefore reduction in the effective flame speed leads to 
reduction in the number of elements. 

The grid evolution with time and adaption level, is shown in 
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Time 

Adp 

Elements 

N 

Apos . err 

Com.grd 

a-used 

a-adp 

Rem 

Ref 

Stf 
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1 

0 

0 

0 

10 
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2 

4 

1 

4 
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2 

94 

94 
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1 
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0 

0 
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2 

13 
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0 
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Table. 5.2. Flame propagation results for 
h version based a strategy (V = 0.0) 
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Figs. 5 . 26a-b . The grid movement is observed to be uniform with 
respect to time during the equilibrium phase. For the numerical 
computations based on the a- strategy; large time step values 
have been taken as compared to the time steps used for the results 
of Chapter 4. Due to this, during the initial time steps, the 
prediction of regions for refinement is not very accurate. 
However, the situation improves as the marching in time takes 
place. Features such as the growth in the width of flame front 
during the ignition phase, constant width of the front, during 
stabilizied phase etc. are observed, as expected. At every 
location in the burnt region, derefinement process takes place 
immediately after the flame front proceeds away. If derefinement 
is not as effective as refinement, element removal is slow in the 
burnt region, as can be seen from figures at some time instants. 
On the other hand, in the unburnt region, the mesh is always 
coarse due to the coarse initial mesh. 

(b). Error behaviour : 

In order to highlight the improvement in error estimation for 
the a-strategy, the variations of the a posteriori error and the 
exponent a with time have been presented in Fig 5.27 for two 
different uniform meshes having 200 and 400 elements, for the flow 
velocity value of V =5.0. The a posteriori error rapidly 
increases in the beginning due to the highly nonlinear ignition 
process occurring near the x = 1.0 boundary. After these initial 
transients, the a posteriori error undergoes systematic bounded 
oscillations with time as discussed in Chapter 4. The oscillations 
in error can be attributed to the motion of the flame front 
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through the regular array of grid points in the uniform mesh. The 
predicted value of a for the uniform mesh (although it is not used 
here for grid adaption) increases during the ignition phase and 
attains a constant value during the steady propagation of the 
front. In fact the value of a also increases and decreases in 
accordance with the error behavior. Comparing the results of the 
N = 200 and N = 400 meshes, on the refined grid the amplitude of 
the oscillations in the a posteriori error is less and the value 
of a is also less. The smaller value of a can be attributed to 
the reduction in the a posteriori error as well as the asymptotic 
stabilization process in the estimation of the error. 

The a posteriori error, a , solution gradient, and the 
maximum and mean of error indicators on adaptive meshes are 
plotted in Figs. 5.28-5.30 for different flow velocities. The 
value of a used at each time step decreases rapidly in the 
beginning and stabilizes for later times. Within each time step 
the magnitude of a for each adaption fluctuates slightly, due to 
the oscillations in the a posteriori error between the limits 
and 5 2 * T ^ e osc ill at i ons in a posteriori error are caused by the 
following alternate sequence of events : (i) error increases when 

the flame front moves out of an existing optimal grid (ii) it 
decreases back due to the readjustment in the grid to the new 
location of the front after the adaption. The fluctuations in the 
a posteriori error are also reflected in other measures of error 
such as the maximum and the mean error indicators. 

The norm of the solution gradient increases during the 
ignition phase and attains a perfectly constant value during the 
equilibrium phase of the front movement. In the test problem, the 
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norm of the gradient error was plotted because of the availability 
of the exact solution. For the flame propagation study, however, 
the norm of the gradient itself has been plotted since the norm of 
gradient error can not be estimated. The constancy of the 
gradient for later times is adequate proof for the efficiency of 
the a scheme and the suitability of the a posteriori error 
expression proposed in this Chapter. 

The error behaviour is more or less same for all the flow 
velocities except for the following differences. For V = 5.0 
situation, spikes are seen occasionally in the error plots and in 
the evolution of a. These spikes indicate temporary mismatch 
between the front location and the grid distribution due to the 
large time steps used in the calculation. However, the spikes are 
eliminated almost immediately in the next time step. Another 
effective means to obtain a smoother behaviour is to select a 
suitably small value of the time step At for the time marching. 
For the flow velocity of V = 200.0, the timewise oscillations of 
the error measures and the exponent a have much lower frequency. 
This implies that the flame movement is very slow in the presence 
of convection flow as discussed earlier. 

In the Figs 5.31 - 5.36, the spatial distribution of error 
indicators are shown at various times for the situations of V = 
0.0 and V = 5.0. The features discussed earlier in sections 
5.6.2. are observed here also. For instance, sharp peaks in the 
error indicator correspond to the flame front location, while wide 
peaks represent inadequately refined elements outside the flame 
zone. The shifting of the sharp peak band with time indicates 
movement of flame front. For the flow velocity of V = 5 . 0 , the 
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wide peaks occur frequently at various time instants. This 
observation substantiates the earlier inference made in connection 
with Fig. 5.29 that perhaps a smaller time step may be necessary 
for eliminating the slight mismatch between the front location and 
the grid distribution. 


5.7.3. Computational time : 

The computational time for the h-version based a strategy 
has been presented in tabular form ( table 5.3). In the algorithm, 
a special case V = 0 (in presence of flow ) is considered and the 
efficiency of the adaptive method is investigated. Various 
quantities have been evaluated in a similar way as described in 
Chapter 4. Further, the source vector, the elementary matrices in 
the upwind formulation and the evaluation of a at each adaption 
for a given time may additionally contribute to the CPU time in 
the algorithm. 

Initially the adaption is carried out starting with 10 
uniform elements. Before ignition, the solution obtained by a 
strategy is very oscillatory because the mesh is coarse and larger 
time steps have been chosen. Consequently, the prediction of a is 
inaccura'te and is a large value. Thus, the number of non-linear 
iterations performed is more in comparison to h version algorithm 
of Chapter 4. Due to large value of oc, the adaption does not take 
place at every time step. As also, the less number of elements 
are required to satisfy the error criterion. Thus the CPU time 
requirements are less at time instants (t = 0 . 00001 , t - 0 . 00002 , 
0.00004). It is observed that the bisection iterative method 
which has been employed here for the determination of a takes more 
time for the convergence when the number of elements is large. 
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Dimension 
-less time 

N 

Nodes 

Adp 

Non-lin 

iter 

CPU time/ 
time step 

( in seconds) 

0.000010 

15 

16 

6 

12 

0.319785 

0.000020 

20 

21 

2 

6 

0.246980 

0.000040 

31 

32 

2 

6 

0.358685 

0.000060 

40 

41 

1 

4 

0.352592 

0.000100 

71 

72 

2 

10 

1.247347 

0.000150 

135 

136 

2 

16 

3 . 995194 

0.000180 

220 

221 

1 

18 

4.945962 

0.000520 

275 

276 

2 

30 

16.586933 

0.001100 

272 

273 

1 

14 

7.317469 

0.001710 

307 

308 

2 

28 

18.661804 

0.003250 

431 

432 

1 

14 

12.781802 


Table. 5.3. CPU time of the h version based a strategy 
algorithm for the flame propagation (V = o.O) 
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This is true near the ignition as well as in the flame front 
region. This increases CPU time for each non-adaption/adaption at 
each time step. Further, the dynamic setting of removal of 
tolerance in the refinement and derefinement processes may also 
affect the CPU time. Hence, at time instants ( t= 0.00006, t 
=0.0001, t = 0.000150, 0.000180, t = 0.000520, 0.003250) the CPU 
time is very large although the number of elements is less. But in 
the present problem, the evaluation of upwind elementary 
matrices is not reguired at all. Thus, the h version based cc 
strategy will show better performance if the calculation of 
upwinding matrices, selection of proper time step and the 
optimisation of the code are taken care off. 

5.8. Conclusions : 

The a-strategy described in the present Chapter has potential 
for being applied as a general purpose grid adaption scheme. It 
involves a self correcting error estimator whose asymptotic 
stabilization leads to the correct error estimate for the problem 
at hand. Results predicted by the a scheme match excellently with 
the analytical solution of the test problem having feature of 
front movement. For the non-linear flame problem also, the 
predictions of the a-strategy match with the corresponding 
h-version results for the no flow situation. It can be inferred 
from the behaviour of the various error measures that the error 
estimator of the a-strategy approaches the actual error in an 
asymptotic sense as h tend to zero in the front zone. 
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is developed for the use of a posteriori analysis. However, in 
order to validate the algorithm for adaptive refinement of grid, a 
problem with known analytical solution and possessing certain 
characteristic of the original problem, has been utilized as a 
test case. 

inement algorithm for a system of partial 
differential eguations in one space dimension, is developed. The 
method uses the concept of "eqidistribution " of errors for the 
residual as an adaptive criterion. A simple statistical procedure 
has been used for refinement of grids appropriate regions where 
the solution changes very rapidly. The mesh refinement strategy 
utilizes a posteriori bounds for a semi-discrete finite element 
set up. 

In the hp refinement algorithm , h- and p-refinement 
algorithms have been conjointly applied for same classes of 
problems. Due to large number of grid adaptions and also the 
conversions from p- to h- version meshes, the hp method requires 
more additional calculations in comparison to h-refinement before 
ignition. However, after the ignition, h as well as hp methods 
perform well. Besides, the hp methods requires less number of 
elements to meet the specific tolerance criterion. The performance 
of h method is extremely good even at the initial time steps while 
in hp version, this may not be true as the approximation becomes 
increasingly sensitive to the quality of the initial mesh. Keeping 
processing cost, or storage limitations in mind, hp method 
performs well in comparison to h method only in the front region. 

The algorithm, called dynamic gr id strategy, has been 
applied to the flame problem in the presence of stationary mixture 
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and as well as to the test problem. It is an efficient method to 
solve for the grid and the solution. The advantage in this method 
lies in the fact that the grid distribution readjusts dynamically 
to concentrate points in regions of larger variation of the 
solution. The method performs better in case of front propagation 
problems in comparison to h and hp methods and is generally, 
computationally in expensive. 

(3). In addition to the above algorithms, a generalized h 
refinement based a strategy adaptive algorithm has been 
developed in the Chapter 5. The method has been applied to the 
flame propagation problem in presence of flow and for a test 
problem using upwind FE formulation. The method is based on 
eqidistribution of error strategy involving local mesh size and 
the local residue of the FE solution. The method performs well for 
the transient problems, especially front propagation problems. 
This algorithm is applied to the flame problem in a stationary 
mixture and to a test problem considered in Chapter 4. A 
comparative study shows that this novel method performs reasonably 
well in comparison to other algorithms. 

6.2. Suggestions for future work : 

(!)• The grid generation algorithm developed here can be 
applied to generate structured grid from the unstructured grid. A 
part of the algorithm developed here has been extended by Suzuki 
(1991) to generate structured surface grids based on unstructured 
grids . 


( 2 ) . 


The method can be used to generate three dimensional 
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grids as well. Note that the grid spacing and adaptive grid 
generation can be achieved by taking appropriate control functions 
in the grid generation equations .The iterative method to obtain 
the grid by solving generation equations can be relaxed so that 
the £- and 77 -constant lines can be drawn only once. 

(3) - The h, h-p, and grid speed algorithms based on 
aposteriori error estimates can be used for the general reactive 
diffusive type problems in higher space dimensions. Especially, 
the refinement based on statistical method is very much useful in 
the front propagation problems as the method has a capabilities to 
generate dense meshes in the required zones. In higher dimensions, 
the combination of all adaptive strategies and the refinement 
based on the local error indicator expressions can be applied to 
variety of transient problems. 

(4) . For a general class of non-self adjoint linear and 
non-linear problems, the proposed self adjusting posteriori error 
estimate can be used in the adaptive refinement algorithms. The 
method is easily extendable to higher dimensions. 

(5) . The concept of dynamic and static zone analysis can also 
be introduced to the algorithm dealing with the h-version based a 
strategy. It may result in less expensive computation than the 
present work. 

(6) . For general class of one dimensional transient problems, 
the hp version based oc strategy algorithm can also be developed in 
straight forward way. 
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ADDENDUM 5 

On page 173 the text is modified to include the following 
paragraph at the end of section 4.8.1. of Chapter 4. 

On the other hand, if the front speed is not a constant 
with respect to time, the success of the scheme depends on how 
accurately the front speed can be predicted. If this prediction 
is not accurate, it is quite possible that errors may accumulate, 
due to the discrepancy between the velocities of the front and the 
grid; in course of time, the numerical simulation may become 
unstable because of error accumulation. 


